Changeset 8647


Ignore:
Timestamp:
06/16/11 14:29:52 (14 years ago)
Author:
Mathieu Morlighem
Message:

no more noise_dampening

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

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h

    r8608 r8647  
    184184        CmResponseEnum,
    185185        CmResponsesEnum,
    186         CmNoiseDmpEnum,
     186        CmNoiseDmpEnum, //TO BE DELETED (not used anymore)
    187187        ConstantEnum,
    188188        NumControlsEnum,

  • issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp

    r8601 r8647  
    4343                parameters->AddObject(new DoubleParam(EpsCmEnum,iomodel->eps_cm));
    4444                parameters->AddObject(new DoubleParam(MeanVelEnum,iomodel->meanvel));
    45                 parameters->AddObject(new DoubleParam(CmNoiseDmpEnum,iomodel->cm_noisedmp));
    4645
    4746                parameters->AddObject(new BoolParam(CmGradientEnum,iomodel->cm_gradient));

  • issm/trunk/src/c/objects/Elements/Element.h

    r8643 r8647  
    5050                virtual double RheologyBbarAbsGradient(bool process_units,int weight_index)=0;
    5151                virtual double DragCoefficientAbsGradient(bool process_units,int weight_index)=0;
    52                 virtual double RegularizeInversion(void)=0;
    5352                virtual double SurfaceArea(void)=0;
    5453                virtual void   InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum)=0;

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

    r8646 r8647  
    39863986                        tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
    39873987                        tria->GradjDragGradient(gradient,resp);
     3988                        delete tria->matice; delete tria;
     3989                        break;
     3990                case RheologyBbarAbsGradientEnum:
     3991                        if (!IsOnBed()) return;
     3992                        tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
     3993                        tira->GradjBGradient(gradient,resp);
    39883994                        delete tria->matice; delete tria;
    39893995                        break;
    … …  
    67246730        /*Affect value to the reduced matrix */
    67256731        for(i=0;i<24;i++) *(Pe_reduced+i)=Pi[i]-Pright[i];
    6726 }
    6727 /*}}}*/
    6728 /*FUNCTION Penta::RegularizeInversion {{{1*/
    6729 double Penta::RegularizeInversion(void){
    6730 
    6731         double J;
    6732         Tria* tria=NULL;
    6733 
    6734         /*recover some inputs: */
    6735         int  approximation;
    6736         inputs->GetParameterValue(&approximation,ApproximationEnum);
    6737 
    6738         /*If on water, return 0: */
    6739         if(IsOnWater())return 0;
    6740 
    6741         /*Bail out if this element if:
    6742          * -> Not MacAyeal and not on the surface
    6743          * -> MacAyeal (2d model) and not on bed) */
    6744         if ((approximation!=MacAyealApproximationEnum && !IsOnSurface()) || (approximation==MacAyealApproximationEnum && !IsOnBed())){
    6745                 return 0;
    6746         }
    6747         else if (approximation==MacAyealApproximationEnum){
    6748 
    6749                 /*This element should be collapsed into a tria element at its base. Create this tria element,
    6750                  * and compute RegularizeInversion*/
    6751 
    6752                 /*Depth Average B*/
    6753                 this->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum);
    6754 
    6755                 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face).
    6756                 J=tria->RegularizeInversion();
    6757                 delete tria->matice; delete tria;
    6758 
    6759                 /*delete B average*/
    6760                 this->matice->inputs->DeleteInput(RheologyBbarEnum);
    6761 
    6762                 return J;
    6763         }
    6764         else{
    6765 
    6766                 /*Depth Average B and put it in inputs*/
    6767                 Penta* penta=GetBasalElement();
    6768                 penta->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum);
    6769                 Input* B_input=penta->matice->inputs->GetInput(RheologyBbarEnum);
    6770                 Input* B_copy=(Input*)B_input->copy();
    6771                 this->matice->inputs->AddInput((Input*)B_copy);
    6772 
    6773                 tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face).
    6774                 J=tria->RegularizeInversion();
    6775                 delete tria->matice; delete tria;
    6776 
    6777                 /*delete B average*/
    6778                 this->matice->inputs->DeleteInput(RheologyBbarEnum);
    6779                 penta->matice->inputs->DeleteInput(RheologyBbarEnum);
    6780 
    6781                 return J;
    6782         }
    67836732}
    67846733/*}}}*/

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

    r8643 r8647  
    7979                void   Configure(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters);
    8080                void   SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters);
    81                 double RegularizeInversion(void);
    8281                void   CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs,Vec dg,Vec df);
    8382                void   CreatePVector(Vec pg, Vec pf);

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

    r8646 r8647  
    29302930                        GradjDragGradient(gradient,resp);
    29312931                        break;
     2932                case RheologyBbarAbsGradientEnum:
     2933                        GradjBGradient(gradient,resp);
     2934                        break;
    29322935                default:
    29332936                        _error_("response %s not supported yet",EnumToStringx(responses[resp]));
    29342937        }
     2938}
     2939/*}}}*/
     2940/*FUNCTION Tria::GradjBGradient{{{1*/
     2941void  Tria::GradjBGradient(Vec gradient, int weight_index){
     2942
     2943        int        i,ig;
     2944        int        doflist1[NUMVERTICES];
     2945        double     Jdet,weight;
     2946        double     xyz_list[NUMVERTICES][3];
     2947        double     dh1dh3[NDOF2][NUMVERTICES];
     2948        double     dk[NDOF2];
     2949        double     grade_g[NUMVERTICES]={0.0};
     2950        GaussTria  *gauss=NULL;
     2951
     2952        /*Retrieve all inputs we will be needing: */
     2953        if(IsOnShelf())return;
     2954        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     2955        GetDofList1(&doflist1[0]);
     2956        Input* rheologyb_input=matice->inputs->GetInput(RheologyBbarEnum); _assert_(rheologyb_input);
     2957        Input* weights_input=inputs->GetInput(WeightsEnum);                _assert_(weights_input);
     2958
     2959        /* Start  looping on the number of gaussian points: */
     2960        gauss=new GaussTria(2);
     2961        for (ig=gauss->begin();ig<gauss->end();ig++){
     2962
     2963                gauss->GaussPoint(ig);
     2964
     2965                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     2966                GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss);
     2967                weights_input->GetParameterValue(&weight,gauss,weight_index);
     2968
     2969                /*Build alpha_complement_list: */
     2970                rheologyb_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
     2971
     2972                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
     2973                for (i=0;i<NUMVERTICES;i++){ grade_g[i]+=-weight*Jdet*gauss->weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
     2974        }
     2975        VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);
     2976
     2977        /*Clean up and return*/
     2978        delete gauss;
    29352979}
    29362980/*}}}*/
    … …  
    29422986        int        doflist[NUMVERTICES];
    29432987        double     vx,vy,lambda,mu,thickness,Jdet;
    2944         double     cm_noisedmp,viscosity_complement;
     2988        double     viscosity_complement;
    29452989        double     dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dB[NDOF2];
    29462990        double     xyz_list[NUMVERTICES][3];
    29472991        double     basis[3],epsilon[3];
    29482992        double     dbasis[NDOF2][NUMVERTICES];
    2949         double     grad_g[NUMVERTICES];
    29502993        double     grad[NUMVERTICES]={0.0};
    29512994        GaussTria *gauss = NULL;
    2952 
    2953         /*retrieve some parameters: */
    2954         this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum);
    29552995
    29562996        /* Get node coordinates and dof list: */
    … …  
    29873027
    29883028                /*standard gradient dJ/dki*/
    2989                 for (i=0;i<NUMVERTICES;i++){
    2990                         grad_g[i]=-viscosity_complement*thickness*( (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1])*Jdet*gauss->weight*basis[i];
    2991                 }
    2992                 /*Add regularization term*/
    2993                 for (i=0;i<NUMVERTICES;i++) grad_g[i]-=cm_noisedmp*Jdet*gauss->weight*(dbasis[0][i]*dB[0]+dbasis[1][i]*dB[1]);
    2994                 for(i=0;i<NUMVERTICES;i++) grad[i]+=grad_g[i];
     3029                for (i=0;i<NUMVERTICES;i++) grad[i]+=-viscosity_complement*thickness*(
     3030                                        (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
     3031                                        )*Jdet*gauss->weight*basis[i];
    29953032        }
    29963033
    … …  
    30083045        int        doflist1[NUMVERTICES];
    30093046        double     vx,vy,lambda,mu,alpha_complement,Jdet;
    3010         double     bed,thickness,Neff,drag,cm_noisedmp;
     3047        double     bed,thickness,Neff,drag;
    30113048        double     xyz_list[NUMVERTICES][3];
    30123049        double     dh1dh3[NDOF2][NUMVERTICES];
    … …  
    30193056        GaussTria  *gauss=NULL;
    30203057
     3058        if(IsOnShelf())return;
     3059
    30213060        /*retrive parameters: */
    30223061        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    3023 
    3024         /*retrieve some parameters ands return if iceshelf: */
    3025         this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum);
    3026         if(IsOnShelf())return;
    30273062        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    30283063        GetDofList1(&doflist1[0]);
    … …  
    30623097                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    30633098                for (i=0;i<NUMVERTICES;i++){
    3064 
    3065                         //standard term dJ/dki
    30663099                        grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss->weight*l1l2l3[i];
    3067 
    3068                         //noise dampening d/dki(1/2*(dk/dx)^2)
    3069                         grade_g_gaussian[i]+=-cm_noisedmp*Jdet*gauss->weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
    30703100                }
    30713101               
    … …  
    30923122        double     bed,thickness,Neff;
    30933123        double     lambda,mu,xi,Jdet,vx,vy,vz;
    3094         double     alpha_complement,drag,cm_noisedmp;
     3124        double     alpha_complement,drag;
    30953125        double     surface_normal[3],bed_normal[3];
    30963126        double     xyz_list[NUMVERTICES][3];
    … …  
    31173147        Input* adjointz_input=inputs->GetInput(AdjointzEnum);        _assert_(adjointz_input);
    31183148
    3119         /*retrieve some parameters: */
    3120         this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum);
    3121 
    31223149        /*Get out if shelf*/
    31233150        if(IsOnShelf())return;
    … …  
    31793206                                                -xi    *(2*drag*alpha_complement*(-vx*bed_normal[0]*bed_normal[2]-vy*bed_normal[1]*bed_normal[2]))
    31803207                                                )*Jdet*gauss->weight*l1l2l3[i];
    3181 
    3182                         //Add regularization term
    3183                         grade_g_gaussian[i]+= - cm_noisedmp*Jdet*gauss->weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
    31843208                }
    31853209
    … …  
    32613285        int        i,ig;
    32623286        int        doflist1[NUMVERTICES];
    3263         double     thickness,Jdet,cm_noisedmp;
     3287        double     thickness,Jdet;
    32643288        double     l1l2l3[3];
    32653289        double     dbasis[NDOF2][NUMVERTICES];
    … …  
    32743298
    32753299        /*Retrieve all inputs we will be needing: */
    3276         this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum);
    32773300        Input* adjoint_input=inputs->GetInput(AdjointEnum);     _assert_(adjoint_input);
    32783301        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    … …  
    32933316
    32943317                for(i=0;i<NUMVERTICES;i++) grade_g[i]+=thickness*Dlambda[0]*Jdet*gauss->weight*l1l2l3[i];
    3295 
    3296                 //noise dampening d/dki(1/2*(dk/dx)^2)
    3297                 //for (i=0;i<NUMVERTICES;i++) grade_g[i]-=cm_noisedmp*Jdet*gauss->weight*(dbasis[0][i]*dp[0]+dbasis[1][i]*dp[1]);
    32983318        }
    32993319
    … …  
    33103330        int        i,ig;
    33113331        int        doflist1[NUMVERTICES];
    3312         double     thickness,Jdet,cm_noisedmp;
     3332        double     thickness,Jdet;
    33133333        double     l1l2l3[3];
    33143334        double     dbasis[NDOF2][NUMVERTICES];
    … …  
    33193339
    33203340        /* Get node coordinates and dof list: */
    3321         this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum);
    33223341        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    33233342        GetDofList1(&doflist1[0]);
    … …  
    33423361
    33433362                for(i=0;i<NUMVERTICES;i++) grade_g[i]+=thickness*Dlambda[1]*Jdet*gauss->weight*l1l2l3[i];
    3344 
    3345                 //noise dampening d/dki(1/2*(dk/dx)^2)
    3346                 //for (i=0;i<NUMVERTICES;i++) grade_g[i]-=cm_noisedmp*Jdet*gauss->weight*(dbasis[0][i]*dp[0]+dbasis[1][i]*dp[1]);
    3347         }
    3348 
     3363        }
    33493364        VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);
    33503365
    … …  
    47294744}
    47304745/*}}}*/
    4731 /*FUNCTION Tria::RegularizeInversion {{{1*/
    4732 double Tria::RegularizeInversion(){
    4733 
    4734         /*constants: */
    4735         const int    numdof=NDOF2*NUMVERTICES;
    4736 
    4737         /* Intermediaries */
    4738         int        i,j,ig;
    4739         int        num_controls;
    4740         int       *control_type=NULL;
    4741         double     Jelem = 0;
    4742         double     cm_noisedmp;
    4743         double     Jdet;
    4744         double     xyz_list[NUMVERTICES][3];
    4745         double     dk[NDOF2],dB[NDOF2];
    4746         GaussTria *gauss = NULL;
    4747 
    4748         /*retrieve parameters and inputs*/
    4749         this->parameters->FindParam(&num_controls,NumControlsEnum);
    4750         this->parameters->FindParam(&control_type,NULL,ControlTypeEnum);
    4751         this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum);
    4752 
    4753         /*If on water, return 0: */
    4754         if(IsOnWater()) return 0;
    4755 
    4756         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    4757 
    4758         for(i=0;i<num_controls;i++){
    4759                 /* Start looping on the number of gaussian points: */
    4760                 gauss=new GaussTria(2);
    4761                 for (ig=gauss->begin();ig<gauss->end();ig++){
    4762 
    4763                         gauss->GaussPoint(ig);
    4764 
    4765                         GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    4766 
    4767                         /*Add Tikhonov regularization term to misfit:
    4768                          *
    4769                          * WARNING: for now, the regularization is only taken into account by the gradient
    4770                          * the misfit reflects the response only!
    4771                          *
    4772                          * */
    4773                         if (control_type[i]==DragCoefficientEnum){
    4774                                 Input* drag_input=inputs->GetInput(DragCoefficientEnum);      _assert_(drag_input);
    4775                                 drag_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
    4776                                 //Jelem+=cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss->weight;
    4777                         }
    4778                         else if (control_type[i]==RheologyBbarEnum){
    4779                                 //nothing
    4780                         }
    4781                         else if (control_type[i]==DhDtEnum){
    4782                                 //nothing
    4783                         }
    4784                         else if (control_type[i]==VxEnum){
    4785                                 //nothing
    4786                         }
    4787                         else if (control_type[i]==VyEnum){
    4788                                 //nothing
    4789                         }
    4790                         else{
    4791                                 _error_("unsupported control type: %s",EnumToStringx(control_type[i]));
    4792                         }
    4793                 }
    4794 
    4795                 delete gauss;
    4796         }
    4797 
    4798         /*Clean up and return*/
    4799         xfree((void**)&control_type);
    4800         return Jelem;
    4801 }
    4802 /*}}}*/
    48034746/*FUNCTION Tria::RheologyBbarAbsGradient{{{1*/
    48044747double Tria::RheologyBbarAbsGradient(bool process_units,int weight_index){

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

    r8646 r8647  
    7676                void   Configure(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters);
    7777                void   SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters);
    78                 double RegularizeInversion(void);
    7978                void   CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs,Vec dg,Vec df);
    8079                void   CreatePVector(Vec pg, Vec pf);

  • issm/trunk/src/c/objects/IoModel.cpp

    r8592 r8647  
    181181        IoModelFetchData(&this->eps_cm,iomodel_handle,"eps_cm");
    182182        IoModelFetchData(&this->tolx,iomodel_handle,"tolx");
    183         IoModelFetchData(&this->cm_noisedmp,iomodel_handle,"cm_noisedmp");
    184183        IoModelFetchData(&this->cm_gradient,iomodel_handle,"cm_gradient");
    185184        IoModelFetchData(&this->eps_res,iomodel_handle,"eps_res");
    … …  
    354353        this->tolx=0;
    355354        this->maxiter=NULL;
    356         this->cm_noisedmp=0;
    357355        this->cm_min=NULL;
    358356        this->cm_max=NULL;

  • issm/trunk/src/c/objects/IoModel.h

    r8592 r8647  
    125125                double  stokesreconditioning;
    126126                int     shelf_dampening;
    127                 double  cm_noisedmp;
    128127                double* cm_min;
    129128                double* cm_max;
    130129                int     cm_gradient;;
    131 
    132                 double  cm_noisedampening;
    133130
    134131                /*control methods: */
Note: See TracChangeset for help on using the changeset viewer.