Changeset 5640


Ignore:
Timestamp:
09/01/10 07:59:46 (15 years ago)
Author:
Mathieu Morlighem
Message:

removed some linearizations

Location:
issm/trunk/src/c
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Container/Inputs.cpp

    r5537 r5640  
    7171}
    7272/*}}}*/
     73/*FUNCTION Inputs::GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type){{{1*/
     74void Inputs::GetParameterValue(double* pvalue,GaussTria* gauss, int enum_type){
     75
     76        vector<Object*>::iterator object;
     77        Input* input=NULL;
     78        bool   found=false;
     79
     80        /*Go through inputs and check whether any input with the same name is already in: */
     81        for ( object=objects.begin() ; object < objects.end(); object++ ){
     82
     83                input=(Input*)(*object);
     84                if (input->EnumType()==enum_type){
     85                        found=true;
     86                        break;
     87                }
     88        }
     89
     90        if (!found){
     91                /*we could not find an input with the correct enum type. No defaults values were provided,
     92                 * error out: */
     93                ISSMERROR("could not find input with enum type %i (%s)",enum_type,EnumToString(enum_type));
     94        }
     95
     96        /*Ok, we have an input if we made it here, request the input to return the values: */
     97        input->GetParameterValue(pvalue,gauss);
     98
     99}
     100/*}}}*/
    73101/*FUNCTION Inputs::GetParameterValue(double* pvalue,double* gauss,int enum_type,double defaultvalue){{{1*/
    74102void Inputs::GetParameterValue(double* pvalue,double* gauss, int enum_type,double defaultvalue){

  • issm/trunk/src/c/Container/Inputs.h

    r5537 r5640  
    1616class Input;
    1717class Node;
     18class GaussTria;
    1819
    1920class Inputs: public DataSet{
    … …  
    4344                void GetParameterValue(double* pvalue,int enum_type);
    4445                void GetParameterValue(double* pvalue,double* gauss,int enum_type);
     46                void GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type);
    4547                void GetParameterValue(double* pvalue,double* gauss,int enum_type,double defaultvalue);
    4648                void GetParameterValues(double* values,double* gauss_pointers, int numgauss,int enum_type);

  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r5635 r5640  
    21582158        switch(analysis_type){
    21592159
    2160                 case DiagnosticHorizAnalysisEnum: case DiagnosticVertAnalysisEnum:
     2160                case DiagnosticHorizAnalysisEnum:
    21612161
    21622162                        /*default vx,vy and vz: either observation or 0 */
    21632163                        if(!iomodel->vx){
    2164                                 if (false && iomodel->vx_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx_obs[penta_vertex_ids[i]-1]/iomodel->yts; // false TO BE DELETED (for NR CM)
     2164                                if (iomodel->vx_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx_obs[penta_vertex_ids[i]-1]/iomodel->yts;
    21652165                                else                 for(i=0;i<6;i++)nodeinputs[i]=0;
    21662166                                this->inputs->AddInput(new PentaVertexInput(VxEnum,nodeinputs));
    … …  
    21692169                        }
    21702170                        if(!iomodel->vy){
    2171                                 if (false && iomodel->vy_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy_obs[penta_vertex_ids[i]-1]/iomodel->yts; // false TO BE DELETED (for NR CM)
     2171                                if (iomodel->vy_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy_obs[penta_vertex_ids[i]-1]/iomodel->yts;
    21722172                                else                 for(i=0;i<6;i++)nodeinputs[i]=0;
    21732173                                this->inputs->AddInput(new PentaVertexInput(VyEnum,nodeinputs));
    … …  
    45794579        double L[numdof];
    45804580        double l1l2l3[3];
    4581         double alpha2_list[3];
    4582         double basalfriction_list[3]={0.0};
    45834581        double basalfriction;
    45844582        double epsilon[6];

  • issm/trunk/src/c/objects/Elements/Tria.cpp

    r5639 r5640  
    688688                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    689689
    690                 /*Add Tikhonov regularization term to misfit*/
     690                /*Add Tikhonov regularization term to misfit:
     691                 *
     692                 * WARNING: for now, the regularization is only taken into account by the gradient
     693                 * the misfit reflects the response only!
     694                 *
     695                 * */
    691696                if (control_type==DragCoefficientEnum){
    692                         if (!shelf) return 0; //only shelf?... TO BE CHECKED
    693697                        ISSMASSERT(drag_input);
    694698                        drag_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
    695                         Jelem+=cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss->weight;
     699                        //Jelem+=cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss->weight;
    696700                }
    697701                else if (control_type==RheologyBbarEnum){
    698702                        ISSMASSERT(Bbar_input);
    699703                        Bbar_input->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], gauss);
    700                         Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss->weight;
     704                        //Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss->weight;
    701705                }
    702706                else if (control_type==DhDtEnum){
    703707                        ISSMASSERT(dhdt_input);
    704708                        dhdt_input->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], gauss);
    705                         Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss->weight;
     709                        //Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss->weight;
    706710                }
    707711                else{
    … …  
    10291033        /*nodal functions: */
    10301034        double l1l2l3[3];
    1031         double    alpha2complement_list[numvertices];                          //TO BE DELETED
    1032         double    gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
    10331035
    10341036        /* strain rate: */
    … …  
    10681070        friction=new Friction("2d",inputs,matpar,analysis_type);
    10691071
    1070         /*COMPUT alpha2_list (TO BE DELETED)*/
    1071         for(i=0;i<numvertices;i++){
    1072                 friction->GetAlphaComplement(&alpha2complement_list[i],&gaussgrids[i][0],VxEnum,VyEnum);
    1073         }
    1074 
    10751072        /*Retrieve all inputs we will be needing: */
    10761073        adjointx_input=inputs->GetInput(AdjointxEnum);
    … …  
    10871084
    10881085                /*Build alpha_complement_list: */
    1089                 //if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum); // TO BE UNCOMMENTED
    1090                 //else alpha_complement=0;
    1091                 TriaRef::GetParameterValue(&alpha_complement,&alpha2complement_list[0],gauss); // TO BE DELETED
     1086                if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum);
     1087                else alpha_complement=0;
    10921088       
    10931089                /*Recover alpha_complement and k: */
    … …  
    25942590        switch(analysis_type){
    25952591
    2596                 case DiagnosticHorizAnalysisEnum: case DiagnosticVertAnalysisEnum:
     2592                case DiagnosticHorizAnalysisEnum:
    25972593
    25982594                        /*default vx,vy and vz: either observation or 0 */
    25992595                        if(!iomodel->vx){
    2600                                 if (false && iomodel->vx_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx_obs[tria_vertex_ids[i]-1]/iomodel->yts; //false TO BE DELETED (for NR CM)
     2596                                if (iomodel->vx_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx_obs[tria_vertex_ids[i]-1]/iomodel->yts;
    26012597                                else                 for(i=0;i<3;i++)nodeinputs[i]=0;
    26022598                                this->inputs->AddInput(new TriaVertexInput(VxEnum,nodeinputs));
    … …  
    26052601                        }
    26062602                        if(!iomodel->vy){
    2607                                 if (false && iomodel->vy_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy_obs[tria_vertex_ids[i]-1]/iomodel->yts; //false TO BE DELETED (for NR CM)
     2603                                if (iomodel->vy_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy_obs[tria_vertex_ids[i]-1]/iomodel->yts;
    26082604                                else                 for(i=0;i<3;i++)nodeinputs[i]=0;
    26092605                                this->inputs->AddInput(new TriaVertexInput(VyEnum,nodeinputs));
    … …  
    32563252        Friction *friction = NULL;
    32573253        double    alpha2;
    3258         double    alpha2_list[numvertices];                                       //TO BE DELETED
    3259         double    gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
    32603254
    32613255        double MAXSLOPE=.06; // 6 %
    … …  
    32953289        friction=new Friction("2d",inputs,matpar,analysis_type);
    32963290
    3297         /*COMPUT alpha2_list (TO BE DELETED)*/
    3298         for(i=0;i<numvertices;i++){
    3299                 friction->GetAlpha2(&alpha2_list[i],&gaussgrids[i][0],VxEnum,VyEnum,VzEnum);
    3300         }
    3301 
    33023291        /* Start  looping on the number of gaussian points: */
    33033292        gauss=new GaussTria(2);
    … …  
    33073296
    33083297                /*Friction: */
    3309                 // friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT
    3310                 TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss); // TO BE DELETED
     3298                friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
    33113299
    33123300                // If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case,
    … …  
    33883376        Friction *friction = NULL;
    33893377        double    alpha2;
    3390         double    alpha2_list[numvertices];                                       //TO BE DELETED
    3391         double    gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
    33923378
    33933379        double MAXSLOPE=.06; // 6 %
    … …  
    34263412        friction=new Friction("2d",inputs,matpar,analysis_type);
    34273413
    3428         /*COMPUT alpha2_list (TO BE DELETED)*/
    3429         for(i=0;i<numvertices;i++){
    3430                 friction->GetAlpha2(&alpha2_list[i],&gaussgrids[i][0],VxEnum,VyEnum,VzEnum);
    3431         }
    3432 
    34333414        /* Start  looping on the number of gaussian points: */
    34343415        gauss=new GaussTria(2);
    … …  
    34383419
    34393420                /*Friction: */
    3440                 // friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT
    3441                 TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss); // TO BE DELETED
     3421                friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
    34423422
    34433423                // If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case,
    … …  
    35173497        Friction *friction = NULL;
    35183498        double    alpha2;
    3519         double    alpha2_list[numvertices];                                       //TO BE DELETED
    3520         double    gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
    35213499
    35223500        double MAXSLOPE=.06; // 6 %
    … …  
    35553533        friction=new Friction("2d",inputs,matpar,analysis_type);
    35563534
    3557         /*COMPUT alpha2_list (TO BE DELETED)*/
    3558         for(i=0;i<numvertices;i++){
    3559                 friction->GetAlpha2(&alpha2_list[i],&gaussgrids[i][0],VxEnum,VyEnum,VzEnum);
    3560         }
    3561 
    35623535        /* Start  looping on the number of gaussian points: */
    35633536        gauss=new GaussTria(2);
    … …  
    35673540
    35683541                /*Friction: */
    3569                 // friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT
    3570                 TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss); // TO BE DELETED
     3542                friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT
    35713543
    35723544                // If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case,
    … …  
    56075579        double alpha2,vx,vy;
    56085580        double geothermalflux_value;
    5609         double alpha2_list[numvertices];                                    //TO BE DELETED
    5610         double gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
    5611         double vx_list[numvertices]; //TO BE DELETED
    5612         double vy_list[numvertices]; //TO BE DELETED
    5613         double basalfriction_list[numvertices]; //TO BE DELETED
    56145581
    56155582        /* gaussian points: */
    … …  
    56295596
    56305597        /*inputs: */
    5631         Input* vx_input=NULL; //TO BE DELETED
    5632         Input* vy_input=NULL; //TO BE DELETED
    5633         Input* vz_input=NULL; //TO BE DELETED
     5598        Input* vx_input=NULL;
     5599        Input* vy_input=NULL;
     5600        Input* vz_input=NULL;
    56345601        Input* geothermalflux_input=NULL;
    56355602       
    … …  
    56565623        geothermalflux_input=inputs->GetInput(GeothermalFluxEnum);
    56575624
    5658         /*COMPUT alpha2_list and basalfriction_list (TO BE DELETED)*/
    5659         for(i=0;i<numvertices;i++){
    5660                 friction->GetAlpha2(&alpha2_list[i],&gaussgrids[i][0],VxEnum,VyEnum,VzEnum); //TO BE DELETED
    5661         }
    5662         vx_input->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3); //TO BE DELETED
    5663         vy_input->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3); //TO BE DELETED
    5664         for(i=0;i<numvertices;i++){
    5665                 basalfriction_list[i]=alpha2_list[i]*(pow(vx_list[i],(double)2.0)+pow(vy_list[i],(double)2.0)); //TO BE DELETED
    5666         }
    5667        
    56685625        /* Ice/ocean heat exchange flux on ice shelf base */
    56695626
    … …  
    56845641       
    56855642                /*Friction: */
    5686                 //friction->GetAlpha2(&alpha2,&gauss_coord[0],VxEnum,VyEnum,VzEnum);
    5687                 TriaRef::GetParameterValue(&basalfriction,&basalfriction_list[0],gauss); // TO BE DELETED
     5643                friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
    56885644               
    56895645                /*Calculate scalar parameter*/

  • issm/trunk/src/c/objects/Loads/Friction.cpp

    r5103 r5640  
    116116}
    117117/*}}}*/
     118/*FUNCTION Friction::GetAlpha2{{{1*/
     119void Friction::GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){
     120
     121        /*This routine calculates the basal friction coefficient
     122          alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/
     123
     124        /*diverse: */
     125        double  r,s;
     126        double  drag_p, drag_q;
     127        double  gravity,rho_ice,rho_water;
     128        double  Neff;
     129        double  thickness,bed;
     130        double  vx,vy,vz,vmag;
     131        double  drag_coefficient;
     132        double  alpha2;
     133
     134
     135        /*Recover parameters: */
     136        inputs->GetParameterValue(&drag_p,DragPEnum);
     137        inputs->GetParameterValue(&drag_q,DragQEnum);
     138        inputs->GetParameterValue(&thickness, gauss,ThicknessEnum);
     139        inputs->GetParameterValue(&bed, gauss,BedEnum);
     140        inputs->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum);
     141
     142        /*Get material parameters: */
     143        gravity=matpar->GetG();
     144        rho_ice=matpar->GetRhoIce();
     145        rho_water=matpar->GetRhoWater();
     146
     147        //compute r and q coefficients: */
     148        r=drag_q/drag_p;
     149        s=1./drag_p;
     150
     151        //From bed and thickness, compute effective pressure when drag is viscous:
     152        Neff=gravity*(rho_ice*thickness+rho_water*bed);
     153
     154        /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because
     155          the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour
     156          for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should
     157          replace it by Neff=0 (ie, equival it to an ice shelf)*/
     158        if (Neff<0)Neff=0;
     159
     160        if(strcmp(element_type,"2d")==0){
     161                inputs->GetParameterValue(&vx, gauss,vxenum);
     162                inputs->GetParameterValue(&vy, gauss,vyenum);
     163                vmag=sqrt(pow(vx,2)+pow(vy,2));
     164        }
     165        else if (strcmp(element_type,"3d")==0){
     166                inputs->GetParameterValue(&vx, gauss,vxenum);
     167                inputs->GetParameterValue(&vy, gauss,vyenum);
     168                inputs->GetParameterValue(&vz, gauss,vzenum);
     169                vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2));
     170        }
     171        else ISSMERROR("element_type %s not supported yet",element_type);
     172
     173        alpha2=pow(drag_coefficient,2)*pow(Neff,r)*pow(vmag,(s-1));
     174
     175        /*Assign output pointers:*/
     176        *palpha2=alpha2;
     177}
     178/*}}}*/
    118179/*FUNCTION Friction::GetAlphaComplement {{{1*/
    119180void Friction::GetAlphaComplement(double* palpha_complement, double* gauss,int vxenum,int vyenum){
    … …  
    174235}
    175236/*}}}*/
     237/*FUNCTION Friction::GetAlphaComplement {{{1*/
     238void Friction::GetAlphaComplement(double* palpha_complement, GaussTria* gauss,int vxenum,int vyenum){
     239
     240        /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p.
     241         * FrictionGetAlphaComplement is used in control methods on drag, and it computes:
     242         * alpha_complement= Neff ^r * vel ^s*/
     243
     244        /*diverse: */
     245        int     i;
     246        double  Neff;
     247        double  r,s;
     248        double  vx;
     249        double  vy;
     250        double  vmag;
     251        double  drag_p,drag_q;
     252        double  drag_coefficient;
     253        double  bed,thickness;
     254        double  gravity,rho_ice,rho_water;
     255        double  alpha_complement;
     256
     257        /*Recover parameters: */
     258        inputs->GetParameterValue(&drag_p,DragPEnum);
     259        inputs->GetParameterValue(&drag_q,DragQEnum);
     260        inputs->GetParameterValue(&thickness, gauss,ThicknessEnum);
     261        inputs->GetParameterValue(&bed, gauss,BedEnum);
     262        inputs->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum);
     263
     264        /*Get material parameters: */
     265        gravity=matpar->GetG();
     266        rho_ice=matpar->GetRhoIce();
     267        rho_water=matpar->GetRhoWater();
     268
     269
     270        //compute r and q coefficients: */
     271        r=drag_q/drag_p;
     272        s=1./drag_p;
     273
     274        //From bed and thickness, compute effective pressure when drag is viscous:
     275        Neff=gravity*(rho_ice*thickness+rho_water*bed);
     276
     277        /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because
     278          the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour
     279          for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should
     280          replace it by Neff=0 (ie, equival it to an ice shelf)*/
     281        if (Neff<0)Neff=0;
     282
     283        //We need the velocity magnitude to evaluate the basal stress:
     284        inputs->GetParameterValue(&vx, gauss,vxenum);
     285        inputs->GetParameterValue(&vy, gauss,vyenum);
     286        vmag=sqrt(pow(vx,2)+pow(vy,2));
     287
     288        alpha_complement=pow(Neff,r)*pow(vmag,(s-1));
     289
     290        /*Assign output pointers:*/
     291        *palpha_complement=alpha_complement;
     292
     293}
     294/*}}}*/

  • issm/trunk/src/c/objects/Loads/Friction.h

    r4007 r5640  
    2828                void  Echo(void);
    2929                void  GetAlpha2(double* palpha2, double* gauss,int vxenum,int vyenum,int vzenum);
     30                void  GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum);
    3031                void  GetAlphaComplement(double* alpha_complement, double* gauss,int vxenum,int vyenum);
     32                void  GetAlphaComplement(double* alpha_complement, GaussTria* gauss,int vxenum,int vyenum);
    3133
    3234};
Note: See TracChangeset for help on using the changeset viewer.