Changeset 22265


Ignore:
Timestamp:
11/16/17 14:12:51 (7 years ago)
Author:
Mathieu Morlighem
Message:

CHG: temporary commit, last code with Jerome

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

Legend:

Unmodified
Added
Removed

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

    r22178 r22265  
    9090                                case SurfaceAbsMisfitEnum:
    9191                                        for(i=0;i<numnodes;i++) pe->values[i]+=(surfaceobs-surface)*weight*Jdet*gauss->weight*basis[i];
     92                                        break;
     93                                case OmegaAbsGradientEnum:
     94                                        /*Nothing in P vector*/
     95                                        break;
     96                                case EtaDiffEnum:
     97                                        /*Nothing in P vector*/
    9298                                        break;
    9399                                default:
    … …  
    132138                        /*Nothing, \partial J/\partial k = 0*/
    133139                        break;
     140                case OmegaAbsGradientEnum: GradientJOmegaGradient(element,gradient,control_index); break;
     141                case EtaDiffEnum: GradientJEtaDiff(element,gradient,control_index); break;
    134142                default: _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
    135143        }
    … …  
    193201        /*Intermediaries*/
    194202        int         n=3;
    195         IssmDouble  dlambda[2],ds[2],slopex,slopey,slope,omega,Jdet;
     203        IssmDouble  dlambda[2],ds[2],slopex,slopey,slope,omega,Jdet,velobs;
    196204        IssmDouble *xyz_list= NULL;
    197205
    … …  
    212220        Input* surfaceslopex_input = element->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input);
    213221        Input* surfaceslopey_input = element->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input);
    214         IssmDouble rhog    = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
     222        Input* velobs_input        = element->GetInput(InversionVelObsEnum); _assert_(velobs_input);
    215223
    216224        Gauss* gauss=element->NewGauss(2);
    … …  
    226234                surfaceslopex_input->GetInputValue(&slopex,gauss);
    227235                surfaceslopey_input->GetInputValue(&slopey,gauss);
     236                velobs_input->GetInputValue(&velobs,gauss);
    228237
    229238                slope = sqrt(slopex*slopex + slopey*slopey);
     239                //if(slope<1.e-5) slope = 1.e-5;
    230240
    231241                /*Build gradient vector (actually -dJ/da): */
    232242                for(int i=0;i<numvertices;i++){
    233                         ge[i]+= -Jdet*gauss->weight*basis[i]*pow(rhog,n)*exp(omega)*pow(slope,n-1)*(ds[0]*dlambda[0] + ds[1]*dlambda[1]);
     243                        ge[i]+= - Jdet*gauss->weight*basis[i]*velobs/slope*(ds[0]*dlambda[0] + ds[1]*dlambda[1]);
    234244                        _assert_(!xIsNan<IssmDouble>(ge[i]));
    235245                }
    … …  
    243253        xDelete<int>(vertexpidlist);
    244254        delete gauss;
     255}/*}}}*/
     256void           AdjointBalancethickness2Analysis::GradientJOmegaGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     257
     258        /*Intermediaries*/
     259        IssmDouble Jdet,weight;
     260        IssmDouble dk[3];
     261        IssmDouble *xyz_list= NULL;
     262
     263        /*Fetch number of vertices for this finite element*/
     264        int numvertices = element->GetNumberOfVertices();
     265
     266        /*Initialize some vectors*/
     267        IssmDouble* dbasis        = xNew<IssmDouble>(2*numvertices);
     268        IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
     269        int*        vertexpidlist = xNew<int>(numvertices);
     270
     271        /*Retrieve all inputs we will be needing: */
     272        element->GetVerticesCoordinates(&xyz_list);
     273        element->GradientIndexing(&vertexpidlist[0],control_index);
     274        Input* omega_input = element->GetInput(BalancethicknessOmegaEnum); _assert_(omega_input);
     275        Input* weights_input         = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     276
     277        /* Start  looping on the number of gaussian points: */
     278        Gauss* gauss=element->NewGauss(2);
     279        for(int ig=gauss->begin();ig<gauss->end();ig++){
     280                gauss->GaussPoint(ig);
     281
     282                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     283                element->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss);
     284                weights_input->GetInputValue(&weight,gauss,OmegaAbsGradientEnum);
     285
     286                /*Build alpha_complement_list: */
     287                omega_input->GetInputDerivativeValue(&dk[0],xyz_list,gauss);
     288
     289                /*Build gradient vector (actually -dJ/ddrag): */
     290                for(int i=0;i<numvertices;i++){
     291                        ge[i]+=-weight*Jdet*gauss->weight*2*(dk[0]*dk[0] + dk[1]*dk[1])*(dbasis[0*numvertices+i]*dk[0]+dbasis[1*numvertices+i]*dk[1]);
     292                        _assert_(!xIsNan<IssmDouble>(ge[i]));
     293                }
     294        }
     295        gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
     296
     297        /*Clean up and return*/
     298        xDelete<IssmDouble>(xyz_list);
     299        xDelete<IssmDouble>(dbasis);
     300        xDelete<IssmDouble>(ge);
     301        xDelete<int>(vertexpidlist);
     302        delete gauss;
     303
     304}/*}}}*/
     305void           AdjointBalancethickness2Analysis::GradientJEtaDiff(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     306
     307        /*Intermediaries*/
     308        IssmDouble Jdet,weight;
     309        IssmDouble omega,omega0;
     310        IssmDouble *xyz_list= NULL;
     311
     312        /*Fetch number of vertices for this finite element*/
     313        int numvertices = element->GetNumberOfVertices();
     314
     315        /*Initialize some vectors*/
     316        IssmDouble* basis         = xNew<IssmDouble>(numvertices);
     317        IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
     318        int*        vertexpidlist = xNew<int>(numvertices);
     319
     320        /*Retrieve all inputs we will be needing: */
     321        element->GetVerticesCoordinates(&xyz_list);
     322        element->GradientIndexing(&vertexpidlist[0],control_index);
     323        Input* omega_input = element->GetInput(BalancethicknessOmegaEnum);   _assert_(omega_input);
     324        Input* omega0_input = element->GetInput(BalancethicknessOmega0Enum); _assert_(omega0_input);
     325        Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     326
     327        /* Start  looping on the number of gaussian points: */
     328        Gauss* gauss=element->NewGauss(2);
     329        for(int ig=gauss->begin();ig<gauss->end();ig++){
     330                gauss->GaussPoint(ig);
     331
     332                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     333                element->NodalFunctionsP1(basis,gauss);
     334                weights_input->GetInputValue(&weight,gauss,EtaDiffEnum);
     335
     336                /*Build alpha_complement_list: */
     337                omega_input->GetInputValue(&omega,gauss);
     338                omega0_input->GetInputValue(&omega0,gauss);
     339
     340                /*Build gradient vector (actually -dJ/ddrag): */
     341                for(int i=0;i<numvertices;i++){
     342                        ge[i]+=-weight*Jdet*gauss->weight*(omega - omega0)*basis[i];
     343                        _assert_(!xIsNan<IssmDouble>(ge[i]));
     344                }
     345        }
     346        gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
     347
     348        /*Clean up and return*/
     349        xDelete<IssmDouble>(xyz_list);
     350        xDelete<IssmDouble>(basis);
     351        xDelete<IssmDouble>(ge);
     352        xDelete<int>(vertexpidlist);
     353        delete gauss;
     354
    245355}/*}}}*/
    246356void           AdjointBalancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/

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

    r18930 r22265  
    3030                void           GradientJdHdt(Element* element,Vector<IssmDouble>* gradient,int control_index);
    3131                void           GradientJOmega(Element* element,Vector<IssmDouble>* gradient,int control_index);
     32                void           GradientJOmegaGradient(Element* element,Vector<IssmDouble>* gradient,int control_index);
     33                void           GradientJEtaDiff(Element* element,Vector<IssmDouble>* gradient,int control_index);
    3234                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3335                void           UpdateConstraints(FemModel* femmodel);

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

    r20016 r22265  
    176176        if(control_type!=VxEnum &&
    177177                control_type!=VyEnum &&
     178                control_type!=BalancethicknessSpcthicknessEnum &&
    178179                control_type!=BalancethicknessThickeningRateEnum){
    179180                _error_("Control "<<EnumToStringx(control_type)<<" not supported");
    … …  
    192193        /*Deal with second term*/
    193194        switch(control_type){
     195                case BalancethicknessSpcthicknessEnum:   GradientJDirichlet(element,gradient,control_index); break;
    194196                case BalancethicknessThickeningRateEnum: GradientJDhDt(element,gradient,control_index); break;
    195197                case VxEnum:                             GradientJVx(  element,gradient,control_index); break;
    … …  
    201203        xDelete<int>(responses);
    202204
     205}/*}}}*/
     206void           AdjointBalancethicknessAnalysis::GradientJDirichlet(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     207
     208        /*Fetch number of vertices for this finite element*/
     209        int numvertices = element->GetNumberOfVertices();
     210
     211        /*Initialize some vectors*/
     212        IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
     213        IssmDouble* lambda        = xNew<IssmDouble>(numvertices);
     214        int*        vertexpidlist = xNew<int>(numvertices);
     215
     216        /*Retrieve all inputs we will be needing: */
     217        element->GradientIndexing(&vertexpidlist[0],control_index);
     218        element->GetInputListOnVertices(lambda,AdjointEnum);
     219
     220        BalancethicknessAnalysis* analysis = new BalancethicknessAnalysis();
     221        ElementMatrix* Ke = analysis->CreateKMatrix(element);
     222        delete analysis;
     223
     224        /*Transpose and return Ke*/
     225        Ke->Transpose();
     226        _assert_(Ke->nrows == numvertices && Ke->ncols == numvertices);
     227
     228        for(int i=0;i<numvertices;i++){
     229                for(int j=0;j<numvertices;j++){
     230                        ge[i] += Ke->values[i*Ke->ncols + j] * lambda[j];
     231                }
     232                //ge[i]=-lambda[i];
     233                _assert_(!xIsNan<IssmDouble>(ge[i]));
     234        }
     235        gradient->SetValues(numvertices,vertexpidlist,ge,INS_VAL);
     236
     237        /*Clean up and return*/
     238        xDelete<IssmDouble>(ge);
     239        xDelete<IssmDouble>(lambda);
     240        xDelete<int>(vertexpidlist);
     241        delete Ke;
    203242}/*}}}*/
    204243void           AdjointBalancethicknessAnalysis::GradientJDhDt(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/

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

    r18930 r22265  
    2828                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    2929                void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
     30                void           GradientJDirichlet(Element* element,Vector<IssmDouble>* gradient,int control_index);
    3031                void           GradientJDhDt(Element* element,Vector<IssmDouble>* gradient,int control_index);
    3132                void           GradientJVx(Element* element,Vector<IssmDouble>* gradient,int control_index);

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

    r22264 r22265  
    3535        iomodel->FetchDataToInput(elements,"md.slr.sealevel",SealevelEnum,0);
    3636        iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
     37        iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum);
     38        iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum);
    3739        iomodel->FetchDataToInput(elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
    3840        iomodel->FetchDataToInput(elements,"md.smb.mass_balance",SmbMassBalanceEnum);
    3941        iomodel->FetchDataToInput(elements,"md.balancethickness.thickening_rate",BalancethicknessThickeningRateEnum);
    40         iomodel->FetchDataToInput(elements,"md.balancethickness.omega",BalancethicknessOmegaEnum);
    41         iomodel->FetchDataToInput(elements,"md.balancethickness.slopex",SurfaceSlopeXEnum);
    42         iomodel->FetchDataToInput(elements,"md.balancethickness.slopey",SurfaceSlopeYEnum);
    4342
    4443        /*Update elements: */
    … …  
    7170
    7271        /*Intermediaries */
    73         IssmDouble  yts = 365*24*3600.;
    74         IssmDouble  Jdet,vx,vy,vel;
     72        IssmDouble  Jdet,ds[2],slope,velobs,omega;
    7573        IssmDouble* xyz_list = NULL;
    7674
    … …  
    8482        /*Retrieve all inputs and parameters*/
    8583        element->GetVerticesCoordinates(&xyz_list);
    86         Input* vx_input = element->GetInput(SurfaceSlopeXEnum); _assert_(vx_input);
    87         Input* vy_input = element->GetInput(SurfaceSlopeYEnum); _assert_(vy_input);
    88 
    89         /*make sure are diffusivisty is large enough*/
    90         vel = sqrt(vx*vx+vy*vy);
    91         if(sqrt(vx*vx+vy*vy)==0.){
    92                 vx = 0.1/yts;
    93                 vy = 0.1/yts;
    94         }
    95         else if(vel<0.1/yts){
    96                 vx = vx/vel*0.1;
    97                 vy = vy/vel*0.1;
    98 
    99         }
     84        Input* omega_input         = element->GetInput(BalancethicknessOmegaEnum); _assert_(omega_input);
     85        Input* surfaceslopex_input = element->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input);
     86        Input* surfaceslopey_input = element->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input);
     87        Input* velobs_input        = element->GetInput(InversionVelObsEnum); _assert_(velobs_input);
    10088
    10189        /* Start  looping on the number of gaussian points: */
    … …  
    10593                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    10694                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    107                 vx_input->GetInputValue(&vx,gauss);
    108                 vy_input->GetInputValue(&vy,gauss);
     95                surfaceslopex_input->GetInputValue(&ds[0],gauss);
     96                surfaceslopey_input->GetInputValue(&ds[1],gauss);
     97                velobs_input->GetInputValue(&velobs,gauss);
     98                omega_input->GetInputValue(&omega,gauss);
     99
     100                slope = sqrt(ds[0]*ds[0] + ds[1]*ds[1]);
     101                //if(slope<1.e-5) slope = 1.e-5;
    109102
    110103                for(int i=0;i<numnodes;i++){
    111104                        for(int j=0;j<numnodes;j++){
    112                                 Ke->values[i*numnodes+j] += gauss->weight*Jdet*(vx*dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+i]*dbasis[1*numnodes+j]);
     105                                Ke->values[i*numnodes+j] += velobs/slope*omega*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]);
    113106                        }
    114107                }
    … …  
    124117
    125118        /*Intermediaries */
    126         IssmDouble  dhdt[2],mb[2],ms[2],Jdet;
     119        IssmDouble  dhdt,mb,ms,Jdet;
    127120        IssmDouble* xyz_list = NULL;
    128121
    … …  
    148141                element->NodalFunctions(basis,gauss);
    149142
    150                 ms_input->GetInputDerivativeValue(&ms[0],xyz_list,gauss);
    151 
    152                 ms_input->GetInputDerivativeValue(&ms[0],xyz_list,gauss);
    153                 mb_input->GetInputDerivativeValue(&mb[0],xyz_list,gauss);
    154                 dhdt_input->GetInputDerivativeValue(&dhdt[0],xyz_list,gauss);
     143                ms_input->GetInputValue(&ms,gauss);
     144                mb_input->GetInputValue(&mb,gauss);
     145                dhdt_input->GetInputValue(&dhdt,gauss);
    155146
    156147                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(
    157                                         (ms[0]+ms[1]-mb[0]-mb[1]-dhdt[0]-dhdt[1])*basis[i]
     148                                        (ms-mb-dhdt)*basis[i]
    158149                                        );
    159150        }
    … …  
    166157}/*}}}*/
    167158void           Balancethickness2Analysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    168                 element->GetSolutionFromInputsOneDof(solution,ThicknessEnum);
     159                element->GetSolutionFromInputsOneDof(solution,SurfaceEnum);
    169160}/*}}}*/
    170161void           Balancethickness2Analysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
    … …  
    173164void           Balancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    174165
    175                         element->InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
    176 
     166        /*Intermediaries*/
     167        IssmDouble  ds[2],s,b,D;
     168        IssmDouble* xyz_list = NULL;
     169
     170        //element->InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
     171        element->InputUpdateFromSolutionOneDof(solution,SurfaceEnum);
     172
     173        /*Fetch number of vertices and allocate velocity vectors*/
     174        int numvertices = element->GetNumberOfVertices();
     175        IssmDouble* vel_list = xNew<IssmDouble>(numvertices);
     176        IssmDouble* vx_list  = xNew<IssmDouble>(numvertices);
     177        IssmDouble* vy_list  = xNew<IssmDouble>(numvertices);
     178
     179        /*Retrieve all inputs and parameters*/
     180        element->GetVerticesCoordinates(&xyz_list);
     181        Input* D_input   = element->GetInput(BalancethicknessDiffusionCoefficientEnum);
     182        Input* H_input   = element->GetInput(ThicknessEnum);                            _assert_(H_input);
     183        Input* s_input   = element->GetInput(SurfaceEnum);                              _assert_(s_input);
     184        Input* b_input   = element->GetInput(BaseEnum);                                 _assert_(b_input);
     185
     186        /*Calculate velocities*/
     187        Gauss* gauss=element->NewGauss();
     188        for(int iv=0;iv<numvertices;iv++){
     189                gauss->GaussVertex(iv);
     190
     191                if(D_input){
     192                        D_input->GetInputValue(&D,gauss);
     193                }
     194                else{
     195                        D = 0.;
     196                }
     197                b_input->GetInputValue(&b,gauss);
     198                s_input->GetInputValue(&s,gauss);
     199                s_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss);
     200
     201                vx_list[iv] = -1./(s-b)*D*ds[0];
     202                vy_list[iv] = -1./(s-b)*D*ds[1];
     203                vel_list[iv] = sqrt(pow(vx_list[iv],2) + pow(vy_list[iv],2));
     204        }
     205
     206        /*Add vx and vy as inputs to the tria element: */
     207        element->AddInput(VxEnum,vx_list,P1Enum);
     208        element->AddInput(VyEnum,vy_list,P1Enum);
     209        element->AddInput(VelEnum,vel_list,P1Enum);
     210
     211        /*Free ressources:*/
     212        delete gauss;
     213        xDelete<IssmDouble>(vy_list);
     214        xDelete<IssmDouble>(vx_list);
     215        xDelete<IssmDouble>(vel_list);
     216        xDelete<IssmDouble>(xyz_list);
    177217}/*}}}*/
    178218void           Balancethickness2Analysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
Note: See TracChangeset for help on using the changeset viewer.