Changeset 18390


Ignore:
Timestamp:
08/14/14 14:38:41 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: implementing new cost functions and gradient for balancethickness2

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

Legend:

Unmodified
Added
Removed

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

    r18371 r18390  
    108108                                        /*J = (H^2 - Hobs^2)^2*/
    109109                                        //for(i=0;i<numnodes;i++){
    110 
    111                                         //      /*H^2 - Hobs^2*/
    112110                                        //      NUMxH2 = 2.*dbasis[0*numnodes+i]*dphi[0]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2);
    113111                                        //      NUMyH2 = 2.*dbasis[1*numnodes+i]*dphi[1]*(dphi[0]*dphi[0] + dphi[1]*dphi[1] - hu2);
    114112                                        //      DENH2 = vbarobs2*vbarobs2+1.e-14;
    115 
    116                                         //      /*Ubar-Ubar_obs*/
    117                                         //      NUMxUbar = (vyobsbar*dphi[0]*dphi[1] - vxobsbar*dphi[1]*dphi[1])*vbarobs*dbasis[0*numnodes+i];
    118                                         //      NUMyUbar = (vyobsbar*dphi[0]*dphi[0] - vxobsbar*dphi[0]*dphi[1])*vbarobs*dbasis[1*numnodes+i];
    119                                         //      DENUbar  = pow(dphi[0]*dphi[0] + dphi[1]*dphi[1],3./2.)+1.e-14;
    120 
    121113                                        //      pe->values[i]+=(NUMxH2+NUMyH2)/DENH2 *weight*Jdet*gauss->weight;
    122114                                        //}
    123115                                        /*J = phi^2*/
    124116                                        //for(i=0;i<numnodes;i++) pe->values[i]+= phi*basis[i]*weight*Jdet*gauss->weight; //OK
    125                                         /*J = grad J ^2*/
     117                                        /*J = grad phi ^2*/
    126118                                        //for(i=0;i<numnodes;i++) pe->values[i]+= (dphi[0]*dbasis[0*numnodes+i] + dphi[1]*dbasis[1*numnodes+i])*weight*Jdet*gauss->weight; //OK
    127119                                        /*J = (ubar - nux*uobs)^2*/
     120                                        //for(i=0;i<numnodes;i++){
     121                                        //      NUMxUbar = (vyobsbar*dphi[0]*dphi[1] - vxobsbar*dphi[1]*dphi[1])*vbarobs*dbasis[0*numnodes+i];
     122                                        //      NUMyUbar = (vyobsbar*dphi[0]*dphi[0] - vxobsbar*dphi[0]*dphi[1])*vbarobs*dbasis[1*numnodes+i];
     123                                        //      DENUbar  = pow(dphi[0]*dphi[0] + dphi[1]*dphi[1],3./2.)+1.e-14;
     124                                        //      pe->values[i]+=(NUMxUbar-NUMyUbar)/DENUbar *weight*Jdet*gauss->weight;
     125                                        //}
     126                                        /*J = 1/2 (vbar ^ gard(phi))^2*/
    128127                                        for(i=0;i<numnodes;i++){
    129                                                 NUMxUbar = (vyobsbar*dphi[0]*dphi[1] - vxobsbar*dphi[1]*dphi[1])*vbarobs*dbasis[0*numnodes+i];
    130                                                 NUMyUbar = (vyobsbar*dphi[0]*dphi[0] - vxobsbar*dphi[0]*dphi[1])*vbarobs*dbasis[1*numnodes+i];
    131                                                 DENUbar  = pow(dphi[0]*dphi[0] + dphi[1]*dphi[1],3./2.)+1.e-14;
    132                                                 pe->values[i]+=(NUMxUbar-NUMyUbar)/DENUbar *weight*Jdet*gauss->weight;
     128                                                pe->values[i]+= weight*Jdet*gauss->weight*
     129                                                  (nuy*vyobs*dphi[0] - nux*vxobs*dphi[1])*
     130                                                  (nuy*vyobs*dbasis[0*numnodes+i] - nux*vxobs*dbasis[1*numnodes+i]);
    133131                                        }
     132
    134133                                        break;
    135134                                default:
    … …  
    168167        element->FindParam(&responses,NULL,InversionCostFunctionsEnum);
    169168
    170         /*Check that control_type is supported*/
    171         if(control_type!=BalancethicknessApparentMassbalanceEnum){
    172                 _error_("Control "<<EnumToStringx(control_type)<<" not supported");
    173         }
    174 
    175169        /*Deal with first part (partial derivative a J with respect to k)*/
    176170        for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
    … …  
    182176        switch(control_type){
    183177                case BalancethicknessApparentMassbalanceEnum: GradientJAdot(element,gradient,control_index); break;
     178                case BalancethicknessNuxEnum: GradientJNux(element,gradient,control_index); break;
     179                case BalancethicknessNuyEnum: GradientJNuy(element,gradient,control_index); break;
    184180                default: _error_("control type not supported yet: " << EnumToStringx(control_type));
    185181        }
    … …  
    234230        delete gauss;
    235231}/*}}}*/
     232void AdjointBalancethickness2Analysis::GradientJNux(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     233
     234        /*Intermediaries*/
     235        IssmDouble Jdet,weight;
     236        IssmDouble vxobs,vyobs;
     237        IssmDouble nux,nuy,dphi[2];
     238        IssmDouble *xyz_list= NULL;
     239
     240        /*Fetch number of vertices for this finite element*/
     241        int numvertices = element->GetNumberOfVertices();
     242
     243        /*Initialize some vectors*/
     244        IssmDouble* basis         = xNew<IssmDouble>(numvertices);
     245        IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
     246        int*        vertexpidlist = xNew<int>(numvertices);
     247
     248        /*Retrieve all inputs we will be needing: */
     249        element->GetVerticesCoordinates(&xyz_list);
     250        element->GradientIndexing(&vertexpidlist[0],control_index);
     251        Input* weights_input   = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     252        Input* vxobs_input     = element->GetInput(BalancethicknessVxObsEnum); _assert_(vxobs_input);
     253        Input* vyobs_input     = element->GetInput(BalancethicknessVyObsEnum); _assert_(vyobs_input);
     254        Input* nux_input       = element->GetInput(BalancethicknessNuxEnum);   _assert_(nux_input);
     255        Input* nuy_input       = element->GetInput(BalancethicknessNuyEnum);   _assert_(nuy_input);
     256        Input* potential_input = element->GetInput(PotentialEnum);             _assert_(potential_input);
     257
     258
     259        Gauss* gauss=element->NewGauss(2);
     260        for(int ig=gauss->begin();ig<gauss->end();ig++){
     261                gauss->GaussPoint(ig);
     262
     263                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     264                element->NodalFunctionsP1(basis,gauss);
     265                weights_input->GetInputValue(&weight,gauss,Balancethickness2MisfitEnum);
     266                vxobs_input->GetInputValue(&vxobs,gauss);
     267                vyobs_input->GetInputValue(&vyobs,gauss);
     268                nux_input->GetInputValue(&nux,gauss);
     269                nuy_input->GetInputValue(&nuy,gauss);
     270                potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
     271
     272                /*Build gradient vector (actually -dJ/da): */
     273                for(int i=0;i<numvertices;i++){
     274                        ge[i]+= - weight*Jdet*gauss->weight*(-vxobs)*dphi[1]*(nuy*vyobs*dphi[0] - nux*vxobs*dphi[1])*basis[i];
     275                        _assert_(!xIsNan<IssmDouble>(ge[i]));
     276                }
     277        }
     278        gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
     279
     280        /*Clean up and return*/
     281        xDelete<IssmDouble>(ge);
     282        xDelete<IssmDouble>(xyz_list);
     283        xDelete<IssmDouble>(basis);
     284        xDelete<int>(vertexpidlist);
     285        delete gauss;
     286}/*}}}*/
     287void AdjointBalancethickness2Analysis::GradientJNuy(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     288
     289        /*Intermediaries*/
     290        IssmDouble Jdet,weight;
     291        IssmDouble vxobs,vyobs;
     292        IssmDouble nux,nuy,dphi[2];
     293        IssmDouble *xyz_list= NULL;
     294
     295        /*Fetch number of vertices for this finite element*/
     296        int numvertices = element->GetNumberOfVertices();
     297
     298        /*Initialize some vectors*/
     299        IssmDouble* basis         = xNew<IssmDouble>(numvertices);
     300        IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
     301        int*        vertexpidlist = xNew<int>(numvertices);
     302
     303        /*Retrieve all inputs we will be needing: */
     304        element->GetVerticesCoordinates(&xyz_list);
     305        element->GradientIndexing(&vertexpidlist[0],control_index);
     306        Input* weights_input   = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     307        Input* vxobs_input     = element->GetInput(BalancethicknessVxObsEnum); _assert_(vxobs_input);
     308        Input* vyobs_input     = element->GetInput(BalancethicknessVyObsEnum); _assert_(vyobs_input);
     309        Input* nux_input       = element->GetInput(BalancethicknessNuxEnum);   _assert_(nux_input);
     310        Input* nuy_input       = element->GetInput(BalancethicknessNuyEnum);   _assert_(nuy_input);
     311        Input* potential_input = element->GetInput(PotentialEnum);             _assert_(potential_input);
     312
     313
     314        Gauss* gauss=element->NewGauss(2);
     315        for(int ig=gauss->begin();ig<gauss->end();ig++){
     316                gauss->GaussPoint(ig);
     317
     318                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     319                element->NodalFunctionsP1(basis,gauss);
     320                weights_input->GetInputValue(&weight,gauss,Balancethickness2MisfitEnum);
     321                vxobs_input->GetInputValue(&vxobs,gauss);
     322                vyobs_input->GetInputValue(&vyobs,gauss);
     323                nux_input->GetInputValue(&nux,gauss);
     324                nuy_input->GetInputValue(&nuy,gauss);
     325                potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
     326
     327                /*Build gradient vector (actually -dJ/da): */
     328                for(int i=0;i<numvertices;i++){
     329                        ge[i]+= - weight*Jdet*gauss->weight*(+vyobs)*dphi[0]*(nuy*vyobs*dphi[0] - nux*vxobs*dphi[1])*basis[i];
     330                        _assert_(!xIsNan<IssmDouble>(ge[i]));
     331                }
     332        }
     333        gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
     334
     335        /*Clean up and return*/
     336        xDelete<IssmDouble>(ge);
     337        xDelete<IssmDouble>(xyz_list);
     338        xDelete<IssmDouble>(basis);
     339        xDelete<int>(vertexpidlist);
     340        delete gauss;
     341}/*}}}*/
    236342void AdjointBalancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    237343        element->InputUpdateFromSolutionOneDof(solution,AdjointEnum);

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

    r18061 r18390  
    2929                void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
    3030                void GradientJAdot(Element* element,Vector<IssmDouble>* gradient,int control_index);
     31                void GradientJNux(Element* element,Vector<IssmDouble>* gradient,int control_index);
     32                void GradientJNuy(Element* element,Vector<IssmDouble>* gradient,int control_index);
    3133                void InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3234                void UpdateConstraints(FemModel* femmodel);

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

    r18371 r18390  
    13611361                        //J +=.5*(dpotential[0]*dpotential[0] + dpotential[1]*dpotential[1])*weight*Jdet*gauss->weight;
    13621362                        /*J = (ubar - nux*uobs)^2*/
    1363                         J +=0.5*((vxbarobs - vxbar)*(vxbarobs - vxbar) + (vybarobs - vybar)*(vybarobs - vybar))*weight*Jdet*gauss->weight;
     1363                        //J +=0.5*((vxbarobs - vxbar)*(vxbarobs - vxbar) + (vybarobs - vybar)*(vybarobs - vybar))*weight*Jdet*gauss->weight;
     1364                        /*J = 1/2 (vbar ^ gard(phi))^2*/
     1365                        J += 0.5*(nuy*vyobs*dpotential[0] - nux*vxobs*dpotential[1])*(nuy*vyobs*dpotential[0] - nux*vxobs*dpotential[1])*weight*Jdet*gauss->weight;
    13641366                }
    13651367
Note: See TracChangeset for help on using the changeset viewer.