Changeset 17975


Ignore:
Timestamp:
05/09/14 14:51:29 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: moving misfits to modules

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

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r17926 r17975  
    277277                virtual void   Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index)=0;
    278278                virtual IssmDouble ThicknessAbsMisfit(void)=0;
    279                 virtual IssmDouble SurfaceAbsVelMisfit(void)=0;
    280                 virtual IssmDouble SurfaceRelVelMisfit(void)=0;
    281                 virtual IssmDouble SurfaceLogVelMisfit(void)=0;
    282279                virtual IssmDouble SurfaceLogVxVyMisfit(void)=0;
    283280                virtual IssmDouble SurfaceAverageVelMisfit(void)=0;
    … …  
    287284                virtual IssmDouble BalancethicknessMisfit(void)=0;
    288285                virtual IssmDouble RheologyBbarAbsGradient(void)=0;
    289                 virtual IssmDouble DragCoefficientAbsGradient(void)=0;
    290286                virtual void   ControlInputGetGradient(Vector<IssmDouble>* gradient,int enum_type,int control_index)=0;
    291287                virtual void   ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index)=0;

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r17967 r17975  
    33653365}
    33663366/*}}}*/
    3367 IssmDouble Penta::SurfaceAbsVelMisfit(void){/*{{{*/
    3368 
    3369         int    approximation;
    3370         IssmDouble J;
    3371         Tria*  tria=NULL;
    3372 
    3373         /*retrieve inputs :*/
    3374         inputs->GetInputValue(&approximation,ApproximationEnum);
    3375 
    3376         /*If on water, return 0: */
    3377         if(!IsIceInElement())return 0;
    3378 
    3379         /*Bail out if this element if:
    3380          * -> Non SSA and not on the surface
    3381          * -> SSA (2d model) and not on bed) */
    3382         if ((approximation!=SSAApproximationEnum && !IsOnSurface()) || (approximation==SSAApproximationEnum && !IsOnBase())){
    3383                 return 0;
    3384         }
    3385         else if (approximation==SSAApproximationEnum){
    3386 
    3387                 /*This element should be collapsed into a tria element at its base. Create this tria element,
    3388                  * and compute SurfaceAbsVelMisfit*/
    3389                 tria=(Tria*)SpawnTria(0,1,2);
    3390                 J=tria->SurfaceAbsVelMisfit();
    3391                 delete tria->material; delete tria;
    3392                 return J;
    3393         }
    3394         else{
    3395 
    3396                 tria=(Tria*)SpawnTria(3,4,5);
    3397                 J=tria->SurfaceAbsVelMisfit();
    3398                 delete tria->material; delete tria;
    3399                 return J;
    3400         }
    3401 }
    3402 /*}}}*/
    3403 IssmDouble Penta::SurfaceLogVelMisfit(void){/*{{{*/
    3404 
    3405         int    approximation;
    3406         IssmDouble J;
    3407         Tria*  tria=NULL;
    3408 
    3409         /*retrieve inputs :*/
    3410         inputs->GetInputValue(&approximation,ApproximationEnum);
    3411 
    3412         /*If on water, return 0: */
    3413         if(!IsIceInElement())return 0;
    3414 
    3415         /*Bail out if this element if:
    3416          * -> Non SSA and not on the surface
    3417          * -> SSA (2d model) and not on bed) */
    3418         if ((approximation!=SSAApproximationEnum && !IsOnSurface()) || (approximation==SSAApproximationEnum && !IsOnBase())){
    3419                 return 0;
    3420         }
    3421         else if (approximation==SSAApproximationEnum){
    3422 
    3423                 /*This element should be collapsed into a tria element at its base. Create this tria element,
    3424                  * and compute SurfaceLogVelMisfit*/
    3425                 tria=(Tria*)SpawnTria(0,1,2); //lower face is 0, upper face is 1.
    3426                 J=tria->SurfaceLogVelMisfit();
    3427                 delete tria->material; delete tria;
    3428                 return J;
    3429         }
    3430         else{
    3431 
    3432                 tria=(Tria*)SpawnTria(3,4,5); //lower face is 0, upper face is 1.
    3433                 J=tria->SurfaceLogVelMisfit();
    3434                 delete tria->material; delete tria;
    3435                 return J;
    3436         }
    3437 }
    3438 /*}}}*/
    34393367IssmDouble Penta::SurfaceLogVxVyMisfit(void){/*{{{*/
    34403368
    … …  
    34753403}
    34763404/*}}}*/
    3477 IssmDouble Penta::SurfaceRelVelMisfit(void){/*{{{*/
     3405IssmDouble Penta::ThicknessAbsGradient(void){/*{{{*/
     3406
     3407        _error_("Not implemented yet");
     3408}
     3409/*}}}*/
     3410IssmDouble Penta::ThicknessAbsMisfit(void){/*{{{*/
    34783411
    34793412        int    approximation;
    … …  
    34863419        /*If on water, return 0: */
    34873420        if(!IsIceInElement())return 0;
    3488 
    3489         /*Bail out if this element if:
    3490          * -> Non SSA and not on the surface
    3491          * -> SSA (2d model) and not on bed) */
    3492         if ((approximation!=SSAApproximationEnum && !IsOnSurface()) || (approximation==SSAApproximationEnum && !IsOnBase())){
    3493                 return 0;
    3494         }
    3495         else if (approximation==SSAApproximationEnum){
    3496 
    3497                 /*This element should be collapsed into a tria element at its base. Create this tria element,
    3498                  * and compute SurfaceRelVelMisfit*/
    3499                 tria=(Tria*)SpawnTria(0,1,2);
    3500                 J=tria->SurfaceRelVelMisfit();
    3501                 delete tria->material; delete tria;
    3502                 return J;
    3503         }
    3504         else{
    3505 
    3506                 tria=(Tria*)SpawnTria(3,4,5);
    3507                 J=tria->SurfaceRelVelMisfit();
    3508                 delete tria->material; delete tria;
    3509                 return J;
    3510         }
    3511 }
    3512 /*}}}*/
    3513 IssmDouble Penta::ThicknessAbsGradient(void){/*{{{*/
    3514 
    3515         _error_("Not implemented yet");
    3516 }
    3517 /*}}}*/
    3518 IssmDouble Penta::ThicknessAbsMisfit(void){/*{{{*/
    3519 
    3520         int    approximation;
    3521         IssmDouble J;
    3522         Tria*  tria=NULL;
    3523 
    3524         /*retrieve inputs :*/
    3525         inputs->GetInputValue(&approximation,ApproximationEnum);
    3526 
    3527         /*If on water, return 0: */
    3528         if(!IsIceInElement())return 0;
    35293421        _error_("Not implemented yet");
    35303422
    35313423        tria=(Tria*)SpawnTria(0,1,2);
    35323424        J=tria->ThicknessAbsMisfit();
    3533         delete tria->material; delete tria;
    3534         return J;
    3535 }
    3536 /*}}}*/
    3537 IssmDouble Penta::DragCoefficientAbsGradient(void){/*{{{*/
    3538 
    3539         IssmDouble J;
    3540         Tria*  tria=NULL;
    3541 
    3542         /*If on water, on shelf or not on bed, skip: */
    3543         if(!IsIceInElement()|| IsFloating() || !IsOnBase()) return 0;
    3544 
    3545         tria=(Tria*)SpawnTria(0,1,2); //lower face is 0, upper face is 1
    3546         J=tria->DragCoefficientAbsGradient();
    35473425        delete tria->material; delete tria;
    35483426        return J;

  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r17926 r17975  
    132132                #endif
    133133
    134                 IssmDouble DragCoefficientAbsGradient(void);
    135134                void   GradientIndexing(int* indexing,int control_index);
    136135                void   Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index);
    … …  
    149148                IssmDouble RheologyBbarAbsGradient(void);
    150149                IssmDouble ThicknessAbsMisfit(void);
    151                 IssmDouble SurfaceAbsVelMisfit(void);
    152                 IssmDouble SurfaceRelVelMisfit(void);
    153                 IssmDouble SurfaceLogVelMisfit(void);
    154150                IssmDouble SurfaceLogVxVyMisfit(void);
    155151                IssmDouble SurfaceAverageVelMisfit(void);

  • issm/trunk-jpl/src/c/classes/Elements/Seg.h

    r17972 r17975  
    168168#endif
    169169
    170                 IssmDouble DragCoefficientAbsGradient(void){_error_("not implemented yet");};
    171170                void       GradientIndexing(int* indexing,int control_index);
    172171                void       Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index){_error_("not implemented yet");};
    … …  
    190189                IssmDouble RheologyBbarAbsGradient(void){_error_("not implemented yet");};
    191190                IssmDouble ThicknessAbsMisfit(void){_error_("not implemented yet");};
    192                 IssmDouble SurfaceAbsVelMisfit(void){_error_("not implemented yet");};
    193191                IssmDouble ThicknessAbsGradient(void){_error_("not implemented yet");};
    194192                IssmDouble ThicknessAlongGradient(void){_error_("not implemented yet");};
    195193                IssmDouble ThicknessAcrossGradient(void){_error_("not implemented yet");};
    196194                IssmDouble BalancethicknessMisfit(void){_error_("not implemented yet");};
    197                 IssmDouble SurfaceRelVelMisfit(void){_error_("not implemented yet");};
    198                 IssmDouble SurfaceLogVelMisfit(void){_error_("not implemented yet");};
    199195                IssmDouble SurfaceLogVxVyMisfit(void){_error_("not implemented yet");};
    200196                IssmDouble SurfaceAverageVelMisfit(void){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Elements/Tetra.h

    r17926 r17975  
    194194                IssmDouble RheologyBbarAbsGradient(void){_error_("not implemented yet");};
    195195                IssmDouble ThicknessAbsMisfit(void){_error_("not implemented yet");};
    196                 IssmDouble SurfaceAbsVelMisfit(void){_error_("not implemented yet");};
    197196                IssmDouble ThicknessAbsGradient(void){_error_("not implemented yet");};
    198197                IssmDouble ThicknessAlongGradient(void){_error_("not implemented yet");};
    199198                IssmDouble ThicknessAcrossGradient(void){_error_("not implemented yet");};
    200199                IssmDouble BalancethicknessMisfit(void){_error_("not implemented yet");};
    201                 IssmDouble SurfaceRelVelMisfit(void){_error_("not implemented yet");};
    202                 IssmDouble SurfaceLogVelMisfit(void){_error_("not implemented yet");};
    203200                IssmDouble SurfaceLogVxVyMisfit(void){_error_("not implemented yet");};
    204201                IssmDouble SurfaceAverageVelMisfit(void){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r17972 r17975  
    35713571}
    35723572/*}}}*/
    3573 IssmDouble Tria::SurfaceLogVelMisfit(void){/*{{{*/
    3574 
    3575         IssmDouble Jelem=0.;
    3576         IssmDouble misfit,Jdet;
    3577         IssmDouble epsvel=2.220446049250313e-16;
    3578         IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/
    3579         IssmDouble velocity_mag,obs_velocity_mag;
    3580         IssmDouble xyz_list[NUMVERTICES][3];
    3581         IssmDouble vx,vy,vxobs,vyobs,weight;
    3582         GaussTria *gauss=NULL;
    3583 
    3584         /*If on water, return 0: */
    3585         if(!IsIceInElement())return 0;
    3586 
    3587         /* Get node coordinates and dof list: */
    3588         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3589 
    3590         /*Retrieve all inputs we will be needing: */
    3591         Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
    3592         Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
    3593         Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
    3594         Input* vxobs_input  =inputs->GetInput(InversionVxObsEnum);     _assert_(vxobs_input);
    3595         Input* vyobs_input  =inputs->GetInput(InversionVyObsEnum);     _assert_(vyobs_input);
    3596 
    3597         /* Start  looping on the number of gaussian points: */
    3598         gauss=new GaussTria(4);
    3599         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3600 
    3601                 gauss->GaussPoint(ig);
    3602 
    3603                 /* Get Jacobian determinant: */
    3604                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3605 
    3606                 /*Get all parameters at gaussian point*/
    3607                 weights_input->GetInputValue(&weight,gauss,SurfaceLogVelMisfitEnum);
    3608                 vx_input->GetInputValue(&vx,gauss);
    3609                 vy_input->GetInputValue(&vy,gauss);
    3610                 vxobs_input->GetInputValue(&vxobs,gauss);
    3611                 vyobs_input->GetInputValue(&vyobs,gauss);
    3612 
    3613                 /*Compute SurfaceLogVelMisfit:
    3614                  *        4         [        vel + eps     ] 2
    3615                  * J = 4 \bar{v}^2 | log ( -----------  ) | 
    3616                  *                 [       vel   + eps    ]
    3617                  *                            obs
    3618                  */
    3619                 velocity_mag    =sqrt(pow(vx,   2)+pow(vy,   2))+epsvel;
    3620                 obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel;
    3621                 misfit=4*pow(meanvel,2)*pow(log(velocity_mag/obs_velocity_mag),2);
    3622 
    3623                 /*Add to cost function*/
    3624                 Jelem+=misfit*weight*Jdet*gauss->weight;
    3625         }
    3626 
    3627         /*clean-up and Return: */
    3628         delete gauss;
    3629         return Jelem;
    3630 }
    3631 /*}}}*/
    36323573IssmDouble Tria::SurfaceLogVxVyMisfit(void){/*{{{*/
    36333574
    … …  
    36893630}
    36903631/*}}}*/
    3691 IssmDouble Tria::SurfaceAbsVelMisfit(void){/*{{{*/
    3692 
    3693         IssmDouble Jelem=0;
    3694         IssmDouble misfit,Jdet;
    3695         IssmDouble vx,vy,vxobs,vyobs,weight;
    3696         IssmDouble xyz_list[NUMVERTICES][3];
    3697         GaussTria *gauss=NULL;
    3698 
    3699         /*If on water, return 0: */
    3700         if(!IsIceInElement())return 0;
    3701 
    3702         /* Get node coordinates and dof list: */
    3703         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3704 
    3705         /*Retrieve all inputs we will be needing: */
    3706         Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
    3707         Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
    3708         Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
    3709         Input* vxobs_input  =inputs->GetInput(InversionVxObsEnum);     _assert_(vxobs_input);
    3710         Input* vyobs_input  =inputs->GetInput(InversionVyObsEnum);     _assert_(vyobs_input);
    3711 
    3712         /* Start  looping on the number of gaussian points: */
    3713         gauss=new GaussTria(2);
    3714         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3715 
    3716                 gauss->GaussPoint(ig);
    3717 
    3718                 /* Get Jacobian determinant: */
    3719                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3720 
    3721                 /*Get all parameters at gaussian point*/
    3722                 weights_input->GetInputValue(&weight,gauss,SurfaceAbsVelMisfitEnum);
    3723                 vx_input->GetInputValue(&vx,gauss);
    3724                 vy_input->GetInputValue(&vy,gauss);
    3725                 vxobs_input->GetInputValue(&vxobs,gauss);
    3726                 vyobs_input->GetInputValue(&vyobs,gauss);
    3727 
    3728                 /*Compute SurfaceAbsVelMisfitEnum:
    3729                  *
    3730                  *      1  [           2              2 ]
    3731                  * J = --- | (u - u   )  +  (v - v   )  |
    3732                  *      2  [       obs            obs   ]
    3733                  *
    3734                  */
    3735                 misfit=0.5*( pow(vx-vxobs,2) + pow(vy-vyobs,2) );
    3736 
    3737                 /*Add to cost function*/
    3738                 Jelem+=misfit*weight*Jdet*gauss->weight;
    3739         }
    3740 
    3741         /*clean up and Return: */
    3742         delete gauss;
    3743         return Jelem;
    3744 }
    3745 /*}}}*/
    3746 IssmDouble Tria::SurfaceRelVelMisfit(void){/*{{{*/
    3747 
    3748         IssmDouble  Jelem=0;
    3749         IssmDouble  scalex=1,scaley=1;
    3750         IssmDouble  misfit,Jdet;
    3751         IssmDouble  epsvel=2.220446049250313e-16;
    3752         IssmDouble  meanvel=3.170979198376458e-05; /*1000 m/yr*/
    3753         IssmDouble  vx,vy,vxobs,vyobs,weight;
    3754         IssmDouble  xyz_list[NUMVERTICES][3];
    3755         GaussTria *gauss=NULL;
    3756 
    3757         /*If on water, return 0: */
    3758         if(!IsIceInElement())return 0;
    3759 
    3760         /* Get node coordinates and dof list: */
    3761         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3762 
    3763         /*Retrieve all inputs we will be needing: */
    3764         Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
    3765         Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
    3766         Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
    3767         Input* vxobs_input  =inputs->GetInput(InversionVxObsEnum);     _assert_(vxobs_input);
    3768         Input* vyobs_input  =inputs->GetInput(InversionVyObsEnum);     _assert_(vyobs_input);
    3769 
    3770         /* Start  looping on the number of gaussian points: */
    3771         gauss=new GaussTria(4);
    3772         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3773 
    3774                 gauss->GaussPoint(ig);
    3775 
    3776                 /* Get Jacobian determinant: */
    3777                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3778 
    3779                 /*Get all parameters at gaussian point*/
    3780                 weights_input->GetInputValue(&weight,gauss,SurfaceRelVelMisfitEnum);
    3781                 vx_input->GetInputValue(&vx,gauss);
    3782                 vy_input->GetInputValue(&vy,gauss);
    3783                 vxobs_input->GetInputValue(&vxobs,gauss);
    3784                 vyobs_input->GetInputValue(&vyobs,gauss);
    3785 
    3786                 /*Compute SurfaceRelVelMisfit:
    3787                  *                       
    3788                  *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
    3789                  * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
    3790                  *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
    3791                  *              obs                        obs                     
    3792                  */
    3793                 scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
    3794                 scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
    3795                 misfit=0.5*(scalex*pow((vx-vxobs),2)+scaley*pow((vy-vyobs),2));
    3796 
    3797                 /*Add to cost function*/
    3798                 Jelem+=misfit*weight*Jdet*gauss->weight;
    3799         }
    3800 
    3801         /*clean up and Return: */
    3802         delete gauss;
    3803         return Jelem;
    3804 }
    3805 /*}}}*/
    38063632IssmDouble Tria::ThicknessAbsGradient(void){/*{{{*/
    38073633
    … …  
    39853811
    39863812        /* clean up and Return: */
    3987         delete gauss;
    3988         return Jelem;
    3989 }
    3990 /*}}}*/
    3991 IssmDouble Tria::DragCoefficientAbsGradient(void){/*{{{*/
    3992 
    3993         /* Intermediaries */
    3994         IssmDouble Jelem = 0;
    3995         IssmDouble weight;
    3996         IssmDouble Jdet;
    3997         IssmDouble xyz_list[NUMVERTICES][3];
    3998         IssmDouble dp[NDOF2];
    3999         GaussTria *gauss = NULL;
    4000 
    4001         /*retrieve parameters and inputs*/
    4002 
    4003         /*If on water, return 0: */
    4004         if(!IsIceInElement()) return 0;
    4005 
    4006         /*Retrieve all inputs we will be needing: */
    4007         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    4008         Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);         _assert_(weights_input);
    4009         Input* drag_input   =inputs->GetInput(FrictionCoefficientEnum); _assert_(drag_input);
    4010 
    4011         /* Start looping on the number of gaussian points: */
    4012         gauss=new GaussTria(2);
    4013         for(int ig=gauss->begin();ig<gauss->end();ig++){
    4014 
    4015                 gauss->GaussPoint(ig);
    4016 
    4017                 /* Get Jacobian determinant: */
    4018                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    4019 
    4020                 /*Get all parameters at gaussian point*/
    4021                 weights_input->GetInputValue(&weight,gauss,DragCoefficientAbsGradientEnum);
    4022                 drag_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
    4023 
    4024                 /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
    4025                 Jelem+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight;
    4026         }
    4027 
    4028         /*Clean up and return*/
    40293813        delete gauss;
    40303814        return Jelem;

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r17926 r17975  
    137137                #endif
    138138
    139                 IssmDouble DragCoefficientAbsGradient(void);
    140139                void       GradientIndexing(int* indexing,int control_index);
    141140                void       Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index);
    … …  
    159158                IssmDouble RheologyBbarAbsGradient(void);
    160159                IssmDouble ThicknessAbsMisfit(void);
    161                 IssmDouble SurfaceAbsVelMisfit(void);
    162160                IssmDouble ThicknessAbsGradient(void);
    163161                IssmDouble ThicknessAlongGradient(void);
    164162                IssmDouble ThicknessAcrossGradient(void);
    165163                IssmDouble BalancethicknessMisfit(void);
    166                 IssmDouble SurfaceRelVelMisfit(void);
    167                 IssmDouble SurfaceLogVelMisfit(void);
    168164                IssmDouble SurfaceLogVxVyMisfit(void);
    169165                IssmDouble SurfaceAverageVelMisfit(void);

  • issm/trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp

    r16314 r17975  
    1010void SurfaceLogVelMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters){
    1111
    12         /*Intermediary*/
    13         int i;
    14         Element* element=NULL;
    15 
    1612        /*output: */
    17         IssmDouble J=0;
     13        IssmDouble J=0.;
    1814        IssmDouble J_sum;
    1915
    2016        /*Compute Misfit: */
    21         for (i=0;i<elements->Size();i++){
    22                 element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
    23                 J+=element->SurfaceLogVelMisfit();
     17        for(int i=0;i<elements->Size();i++){
     18                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
     19                J+=SurfaceLogVelMisfit(element);
    2420        }
    2521
    … …  
    3228        *pJ=J;
    3329}
     30
     31IssmDouble SurfaceLogVelMisfit(Element* element){
     32
     33        int        domaintype,numcomponents;
     34        bool       surfaceintegral;
     35        IssmDouble Jelem=0.;
     36        IssmDouble epsvel=2.220446049250313e-16;
     37        IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/
     38        IssmDouble velocity_mag,obs_velocity_mag;
     39        IssmDouble misfit,Jdet;
     40        IssmDouble vx,vy,vxobs,vyobs,weight;
     41        IssmDouble* xyz_list = NULL;
     42
     43        /*Get basal element*/
     44        if(!element->IsOnSurface()) return 0.;
     45
     46        /*If on water, return 0: */
     47        if(!element->IsIceInElement()) return 0.;
     48
     49        /*Get problem dimension*/
     50        element->FindParam(&domaintype,DomainTypeEnum);
     51        switch(domaintype){
     52                case Domain2DverticalEnum:
     53                        surfaceintegral = true;
     54                        numcomponents   = 1;
     55                        break;
     56                case Domain3DEnum:
     57                        surfaceintegral = true;
     58                        numcomponents   = 2;
     59                        break;
     60                case Domain2DhorizontalEnum:
     61                        surfaceintegral = false;
     62                        numcomponents   = 2;
     63                        break;
     64                default: _error_("not supported yet");
     65        }
     66
     67        /*Spawn surface element*/
     68        Element* topelement = element->SpawnTopElement();
     69
     70        /* Get node coordinates*/
     71        topelement->GetVerticesCoordinates(&xyz_list);
     72
     73        /*Retrieve all inputs we will be needing: */
     74        Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     75        Input* vx_input     =topelement->GetInput(VxEnum);                                 _assert_(vx_input);
     76        Input* vxobs_input  =topelement->GetInput(InversionVxObsEnum);                     _assert_(vxobs_input);
     77        Input* vy_input     = NULL;
     78        Input* vyobs_input  = NULL;
     79        if(numcomponents==2){
     80                vy_input    =topelement->GetInput(VyEnum);              _assert_(vy_input);
     81                vyobs_input =topelement->GetInput(InversionVyObsEnum);  _assert_(vyobs_input);
     82        }
     83
     84        /* Start  looping on the number of gaussian points: */
     85        Gauss* gauss=topelement->NewGauss(4);
     86        for(int ig=gauss->begin();ig<gauss->end();ig++){
     87
     88                gauss->GaussPoint(ig);
     89
     90                /* Get Jacobian determinant: */
     91                topelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     92
     93                /*Get all parameters at gaussian point*/
     94                weights_input->GetInputValue(&weight,gauss,SurfaceLogVelMisfitEnum);
     95                vx_input->GetInputValue(&vx,gauss);
     96                vxobs_input->GetInputValue(&vxobs,gauss);
     97                if(numcomponents==2){
     98                        vy_input->GetInputValue(&vy,gauss);
     99                        vyobs_input->GetInputValue(&vyobs,gauss);
     100                }
     101
     102                /*Compute SurfaceLogVelMisfit:
     103                 *                 [        vel + eps     ] 2
     104                 * J = 4 \bar{v}^2 | log ( -----------  ) | 
     105                 *                 [       vel   + eps    ]
     106                 *                            obs
     107                 */
     108                if(numcomponents==1){
     109                        velocity_mag    =fabs(vx)+epsvel;
     110                        obs_velocity_mag=fabs(vxobs)+epsvel;
     111                }
     112                else{
     113                        velocity_mag    =sqrt(vx*vx+vy*vy)+epsvel;
     114                        obs_velocity_mag=sqrt(vxobs*vxobs+vyobs*vyobs)+epsvel;
     115                }
     116
     117                misfit=4*pow(meanvel,2)*pow(log(velocity_mag/obs_velocity_mag),2);
     118
     119                /*Add to cost function*/
     120                Jelem+=misfit*weight*Jdet*gauss->weight;
     121        }
     122
     123        /*clean up and Return: */
     124        if(domaintype!=Domain2DhorizontalEnum){topelement->DeleteMaterials(); delete topelement;};
     125        xDelete<IssmDouble>(xyz_list);
     126        delete gauss;
     127        return Jelem;
     128}

  • issm/trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.h

    r16314 r17975  
    1010/* local prototypes: */
    1111void SurfaceLogVelMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters);
     12IssmDouble SurfaceLogVelMisfit(Element* element);
    1213
    1314#endif

  • issm/trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp

    r16314 r17975  
    1010void SurfaceRelVelMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters){
    1111
    12         /*Intermediary*/
    13         int i;
    14         Element* element=NULL;
    15 
    1612        /*output: */
    17         IssmDouble J=0;
     13        IssmDouble J=0.;
    1814        IssmDouble J_sum;
    1915
    2016        /*Compute Misfit: */
    21         for (i=0;i<elements->Size();i++){
    22                 element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
    23                 J+=element->SurfaceRelVelMisfit();
     17        for(int i=0;i<elements->Size();i++){
     18                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
     19                J+=SurfaceRelVelMisfit(element);
    2420        }
    2521
    … …  
    3228        *pJ=J;
    3329}
     30
     31IssmDouble SurfaceRelVelMisfit(Element* element){
     32
     33        int        domaintype,numcomponents;
     34        bool       surfaceintegral;
     35        IssmDouble Jelem=0.;
     36        IssmDouble misfit,Jdet,scalex,scaley;
     37        IssmDouble epsvel=2.220446049250313e-16;
     38        IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/
     39        IssmDouble vx,vy,vxobs,vyobs,weight;
     40        IssmDouble* xyz_list = NULL;
     41
     42        /*Get basal element*/
     43        if(!element->IsOnSurface()) return 0.;
     44
     45        /*If on water, return 0: */
     46        if(!element->IsIceInElement()) return 0.;
     47
     48        /*Get problem dimension*/
     49        element->FindParam(&domaintype,DomainTypeEnum);
     50        switch(domaintype){
     51                case Domain2DverticalEnum:
     52                        surfaceintegral = true;
     53                        numcomponents   = 1;
     54                        break;
     55                case Domain3DEnum:
     56                        surfaceintegral = true;
     57                        numcomponents   = 2;
     58                        break;
     59                case Domain2DhorizontalEnum:
     60                        surfaceintegral = false;
     61                        numcomponents   = 2;
     62                        break;
     63                default: _error_("not supported yet");
     64        }
     65
     66        /*Spawn surface element*/
     67        Element* topelement = element->SpawnTopElement();
     68
     69        /* Get node coordinates*/
     70        topelement->GetVerticesCoordinates(&xyz_list);
     71
     72        /*Retrieve all inputs we will be needing: */
     73        Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     74        Input* vx_input     =topelement->GetInput(VxEnum);                                 _assert_(vx_input);
     75        Input* vxobs_input  =topelement->GetInput(InversionVxObsEnum);                     _assert_(vxobs_input);
     76        Input* vy_input     = NULL;
     77        Input* vyobs_input  = NULL;
     78        if(numcomponents==2){
     79                vy_input    =topelement->GetInput(VyEnum);              _assert_(vy_input);
     80                vyobs_input =topelement->GetInput(InversionVyObsEnum);  _assert_(vyobs_input);
     81        }
     82
     83        /* Start  looping on the number of gaussian points: */
     84        Gauss* gauss=topelement->NewGauss(4);
     85        for(int ig=gauss->begin();ig<gauss->end();ig++){
     86
     87                gauss->GaussPoint(ig);
     88
     89                /* Get Jacobian determinant: */
     90                topelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     91
     92                /*Get all parameters at gaussian point*/
     93                weights_input->GetInputValue(&weight,gauss,SurfaceRelVelMisfitEnum);
     94                vx_input->GetInputValue(&vx,gauss);
     95                vxobs_input->GetInputValue(&vxobs,gauss);
     96                if(numcomponents==2){
     97                        vy_input->GetInputValue(&vy,gauss);
     98                        vyobs_input->GetInputValue(&vyobs,gauss);
     99                }
     100
     101                /*Compute SurfaceRelVelMisfit:
     102                 *                       
     103                 *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
     104                 * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
     105                 *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
     106                 *              obs                        obs                     
     107                 */
     108
     109                if(numcomponents==2){
     110                        scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
     111                        scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
     112                        misfit=0.5*(scalex*pow((vx-vxobs),2)+scaley*pow((vy-vyobs),2));
     113                }
     114                else{
     115                        scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
     116                        misfit=0.5*(scalex*pow((vx-vxobs),2));
     117                }
     118
     119                /*Add to cost function*/
     120                Jelem+=misfit*weight*Jdet*gauss->weight;
     121        }
     122
     123        /*clean up and Return: */
     124        if(domaintype!=Domain2DhorizontalEnum){topelement->DeleteMaterials(); delete topelement;};
     125        xDelete<IssmDouble>(xyz_list);
     126        delete gauss;
     127        return Jelem;
     128}

  • issm/trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.h

    r16314 r17975  
    1010/* local prototypes: */
    1111void SurfaceRelVelMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters);
     12IssmDouble SurfaceRelVelMisfit(Element* element);
    1213
    1314#endif
Note: See TracChangeset for help on using the changeset viewer.