Changeset 8643


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

Better way to deal with gradients of models with respect to parameters

Location:
issm/trunk/src/c/objects/Elements
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r8608 r8643  
    4141                virtual void   GetParameterValue(double* pvalue,Node* node,int enumtype)=0;
    4242                virtual void   Gradj(Vec gradient,int control_type)=0;
    43                 virtual void   GradjDrag(Vec gradient)=0;
    44                 virtual void   GradjB(Vec gradient)=0;
    4543                virtual double ThicknessAbsMisfit(bool process_units  ,int weight_index)=0;
    4644                virtual double SurfaceAbsVelMisfit(bool process_units ,int weight_index)=0;

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

    r8611 r8643  
    38773877void  Penta::Gradj(Vec gradient,int control_type){
    38783878
     3879        int              i,approximation;
     3880        Tria*            tria=NULL;
     3881
    38793882        /*If on water, skip grad (=0): */
    38803883        if(IsOnWater())return;
    38813884
    3882         if (control_type==DragCoefficientEnum){
    3883                 GradjDrag(gradient);
    3884         }
    3885         else if (control_type==RheologyBbarEnum){
    3886                 GradjB(gradient);
    3887         }
    3888         else _error_("control type %s not supported yet: ",EnumToStringx(control_type));
    3889 }
    3890 /*}}}*/
    3891 /*FUNCTION Penta::GradjB {{{1*/
    3892 void  Penta::GradjB(Vec gradient){
    3893 
    3894         int              i;
    3895         int              approximation;
    3896         Tria*            tria           =NULL;
    3897         TriaVertexInput* triavertexinput=NULL;
    3898 
    3899         /*If on water, skip: */
    3900         if(IsOnWater())return;
    3901         inputs->GetParameterValue(&approximation,ApproximationEnum);
    3902 
    3903         if (approximation==MacAyealApproximationEnum){
    3904                 /*Bail out element if MacAyeal (2d) and not on bed*/
    3905                 if (!IsOnBed()) return;
    3906 
    3907                 /*This element should be collapsed into a tria element at its base. Create this tria element,
    3908                  * and compute gardj*/
    3909 
    3910                 /*Depth Average B*/
    3911                 this->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum);
    3912 
    3913                 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face).
    3914                 tria->GradjB(gradient);
    3915                 delete tria->matice; delete tria;
    3916 
    3917                 /*delete B average*/
    3918                 this->matice->inputs->DeleteInput(RheologyBbarEnum);
    3919         }
    3920         else{
    3921                 /*Gradient is computed on bed only (Bbar)*/
    3922                 if (!IsOnBed()) return;
    3923 
    3924                 /*Depth Average B*/
    3925                 this->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum);
    3926 
    3927                 /*B is a 2d field, use MacAyeal(2d) gradient even if it is Stokes or Pattyn*/
    3928                 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face).
    3929                 tria->GradjB(gradient);
    3930                 delete tria->matice; delete tria;
    3931 
    3932                 /*delete B average*/
    3933                 this->matice->inputs->DeleteInput(RheologyBbarEnum);
    3934         }
    3935 }
    3936 /*}}}*/
    3937 /*FUNCTION Penta::GradjDrag {{{1*/
    3938 void  Penta::GradjDrag(Vec gradient){
    3939 
    3940         int              i,approximation;
    3941         double           temp_gradient[6]={0,0,0,0,0,0};
    3942         Tria*            tria=NULL;
    3943         TriaVertexInput* triavertexinput=NULL;
    3944 
    3945         /*retrieve inputs :*/
    3946         inputs->GetParameterValue(&approximation,ApproximationEnum);
    3947 
    3948         /*If on water, on shelf or not on bed, skip: */
    3949         if(IsOnWater()|| IsOnShelf() || !IsOnBed())return;
    3950 
    3951         if (approximation==MacAyealApproximationEnum || approximation==PattynApproximationEnum){
    3952                 /*MacAyeal or Pattyn*/
    3953                 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
    3954                 tria->GradjDrag(gradient);
    3955                 delete tria->matice; delete tria;
    3956         }
    3957         else if (approximation==StokesApproximationEnum){
    3958                 /*Stokes*/
    3959                 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
    3960                 tria->GradjDragStokes(gradient);
    3961                 delete tria->matice; delete tria;
    3962         }
    3963         else if (approximation==NoneApproximationEnum){
    3964                 return;
    3965         }
    3966         else _error_("approximation %s not supported yet",EnumToStringx(approximation));
     3885        switch(control_type){
     3886
     3887                case DragCoefficientEnum:
     3888
     3889                        inputs->GetParameterValue(&approximation,ApproximationEnum);
     3890                        if(IsOnShelf() || !IsOnBed()) return;
     3891
     3892                        switch(approximation){
     3893                                case MacAyealApproximationEnum:
     3894                                        tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
     3895                                        tria->GradjDragMacAyeal(gradient);
     3896                                        delete tria->matice; delete tria;
     3897                                        break;
     3898                                case PattynApproximationEnum:
     3899                                        tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
     3900                                        tria->GradjDragMacAyeal(gradient);
     3901                                        delete tria->matice; delete tria;
     3902                                        //GradjDragPattyn(gradient);
     3903                                        break;
     3904                                case StokesApproximationEnum:
     3905                                        tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
     3906                                        tria->GradjDragStokes(gradient);
     3907                                        delete tria->matice; delete tria;
     3908                                        //GradjDragStokes(gradient);
     3909                                        break;
     3910                                case NoneApproximationEnum:
     3911                                        /*Gradient is 0*/
     3912                                        break;
     3913                                default:
     3914                                        _error_("approximation %s not supported yet",EnumToStringx(approximation));
     3915                        }
     3916                        break;
     3917
     3918                case RheologyBbarEnum:
     3919                        inputs->GetParameterValue(&approximation,ApproximationEnum);
     3920                        switch(approximation){
     3921
     3922                                case MacAyealApproximationEnum:
     3923
     3924                                        /*This element should be collapsed into a tria element at its base*/
     3925                                        if (!IsOnBed()) return;
     3926
     3927                                        /*Depth Average B*/
     3928                                        this->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum);
     3929
     3930                                        tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face).
     3931                                        tria->GradjBMacAyeal(gradient);
     3932                                        delete tria->matice; delete tria;
     3933
     3934                                        /*delete B average*/
     3935                                        this->matice->inputs->DeleteInput(RheologyBbarEnum);
     3936                                        break;
     3937                                case PattynApproximationEnum: case StokesApproximationEnum:
     3938
     3939                                        /*Gradient is computed on bed only (Bbar)*/
     3940                                        if (!IsOnBed()) return;
     3941
     3942                                        /*Depth Average B*/
     3943                                        this->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum);
     3944
     3945                                        /*B is a 2d field, use MacAyeal(2d) gradient even if it is Stokes or Pattyn*/
     3946                                        tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face).
     3947                                        tria->GradjBMacAyeal(gradient);
     3948                                        delete tria->matice; delete tria;
     3949
     3950                                        /*delete B average*/
     3951                                        this->matice->inputs->DeleteInput(RheologyBbarEnum);
     3952                                        break;
     3953
     3954                                case NoneApproximationEnum:
     3955                                        /*Gradient is 0*/
     3956                                        break;
     3957                                default:
     3958                                        _error_("approximation %s not supported yet",EnumToStringx(approximation));
     3959                        }
     3960                        break;
     3961
     3962                default:
     3963                        _error_("control type %s not supported yet: ",EnumToStringx(control_type));
     3964        }
    39673965}
    39683966/*}}}*/

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

    r8608 r8643  
    8989                void   GetVectorFromInputs(Vec vector,int NameEnum);
    9090                void   Gradj(Vec gradient,int control_type);
    91                 void   GradjB(Vec gradient);
    92                 void   GradjDrag(Vec gradient);
    9391                int    Sid();
    9492                void   InputControlUpdate(double scalar,bool save_parameter);

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

    r8611 r8643  
    28822882        switch(control_type){
    28832883                case DragCoefficientEnum:
    2884                         GradjDrag(gradient);
     2884                        GradjDragMacAyeal(gradient);
    28852885                        break;
    28862886                case RheologyBbarEnum:
    2887                         GradjB(gradient);
     2887                        GradjBMacAyeal(gradient);
    28882888                        break;
    28892889                case DhDtEnum:
    2890                         GradjDhDt(gradient);
     2890                        GradjDhDtBalancedthickness(gradient);
    28912891                        break;
    28922892                case VxEnum:
    2893                         GradjVx(gradient);
     2893                        GradjVxBalancedthickness(gradient);
    28942894                        break;
    28952895                case VyEnum:
    2896                         GradjVy(gradient);
     2896                        GradjVyBalancedthickness(gradient);
    28972897                        break;
    28982898                default:
    … …  
    29012901}
    29022902/*}}}*/
    2903 /*FUNCTION Tria::GradjB{{{1*/
    2904 void  Tria::GradjB(Vec gradient){
     2903/*FUNCTION Tria::GradjBMacAyeal{{{1*/
     2904void  Tria::GradjBMacAyeal(Vec gradient){
    29052905
    29062906        /*Intermediaries*/
    … …  
    29672967}
    29682968/*}}}*/
    2969 /*FUNCTION Tria::GradjDrag {{{1*/
    2970 void  Tria::GradjDrag(Vec gradient){
     2969/*FUNCTION Tria::GradjDragMacAyeal {{{1*/
     2970void  Tria::GradjDragMacAyeal(Vec gradient){
    29712971
    29722972        int        i,ig;
    … …  
    31603160}
    31613161/*}}}*/
    3162 /*FUNCTION Tria::GradjDhDt{{{1*/
    3163 void  Tria::GradjDhDt(Vec gradient){
     3162/*FUNCTION Tria::GradjDhDtBalancedthickness{{{1*/
     3163void  Tria::GradjDhDtBalancedthickness(Vec gradient){
    31643164
    31653165        /*Intermediaries*/
    … …  
    31773177}
    31783178/*}}}*/
    3179 /*FUNCTION Tria::GradjVx{{{1*/
    3180 void  Tria::GradjVx(Vec gradient){
     3179/*FUNCTION Tria::GradjVxBalancedthickness{{{1*/
     3180void  Tria::GradjVxBalancedthickness(Vec gradient){
    31813181
    31823182        /*Intermediaries*/
    … …  
    32263226}
    32273227/*}}}*/
    3228 /*FUNCTION Tria::GradjVy{{{1*/
    3229 void  Tria::GradjVy(Vec gradient){
     3228/*FUNCTION Tria::GradjVyBalancedthickness{{{1*/
     3229void  Tria::GradjVyBalancedthickness(Vec gradient){
    32303230
    32313231        /*Intermediaries*/

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

    r8608 r8643  
    9090                void   GetVectorFromInputs(Vec vector,int NameEnum);
    9191                void   Gradj(Vec gradient,int control_type);
    92                 void   GradjB(Vec gradient);
    93                 void   GradjDrag(Vec gradient);
    94                 void   GradjDhDt(Vec gradient);
    95                 void   GradjVx(Vec gradient);
    96                 void   GradjVy(Vec gradient);
     92                void   GradjBMacAyeal(Vec gradient);
     93                void   GradjDragMacAyeal(Vec gradient);
     94                void   GradjDragStokes(Vec gradient);
     95                void   GradjDhDtBalancedthickness(Vec gradient);
     96                void   GradjVxBalancedthickness(Vec gradient);
     97                void   GradjVyBalancedthickness(Vec gradient);
    9798                void   InputControlUpdate(double scalar,bool save_parameter);
    9899                void   InputArtificialNoise(int enum_type,double min, double max);
    … …  
    187188                void      GetSolutionFromInputsHydrology(Vec solution);
    188189                void    GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input);
    189                 void      GradjDragStokes(Vec gradient);
    190190                void      InputUpdateFromSolutionAdjointBalancethickness( double* solution);
    191191                void      InputUpdateFromSolutionAdjointHoriz( double* solution);
Note: See TracChangeset for help on using the changeset viewer.