Changeset 17976


Ignore:
Timestamp:
05/09/14 15:37:18 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: pushing misfits away from elements

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

Legend:

Unmodified
Added
Removed

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

    r17975 r17976  
    276276
    277277                virtual void   Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index)=0;
    278                 virtual IssmDouble ThicknessAbsMisfit(void)=0;
    279                 virtual IssmDouble SurfaceLogVxVyMisfit(void)=0;
    280                 virtual IssmDouble SurfaceAverageVelMisfit(void)=0;
    281                 virtual IssmDouble ThicknessAbsGradient(void)=0;
    282                 virtual IssmDouble ThicknessAlongGradient(void)=0;
    283                 virtual IssmDouble ThicknessAcrossGradient(void)=0;
    284                 virtual IssmDouble BalancethicknessMisfit(void)=0;
    285                 virtual IssmDouble RheologyBbarAbsGradient(void)=0;
    286278                virtual void   ControlInputGetGradient(Vector<IssmDouble>* gradient,int enum_type,int control_index)=0;
    287279                virtual void   ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index)=0;

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

    r17975 r17976  
    33293329}
    33303330/*}}}*/
    3331 IssmDouble Penta::SurfaceAverageVelMisfit(void){/*{{{*/
    3332 
    3333         int    approximation;
    3334         IssmDouble J;
    3335         Tria*  tria=NULL;
    3336 
    3337         /*retrieve inputs :*/
    3338         inputs->GetInputValue(&approximation,ApproximationEnum);
    3339 
    3340         /*If on water, return 0: */
    3341         if(!IsIceInElement())return 0;
    3342 
    3343         /*Bail out if this element if:
    3344          * -> Non SSA and not on the surface
    3345          * -> SSA (2d model) and not on bed) */
    3346         if ((approximation!=SSAApproximationEnum && !IsOnSurface()) || (approximation==SSAApproximationEnum && !IsOnBase())){
    3347                 return 0;
    3348         }
    3349         else if (approximation==SSAApproximationEnum){
    3350 
    3351                 /*This element should be collapsed into a tria element at its base. Create this tria element,
    3352                  * and compute SurfaceAverageVelMisfit*/
    3353                 tria=(Tria*)SpawnTria(0,1,2);
    3354                 J=tria->SurfaceAverageVelMisfit();
    3355                 delete tria->material; delete tria;
    3356                 return J;
    3357         }
    3358         else{
    3359 
    3360                 tria=(Tria*)SpawnTria(3,4,5);
    3361                 J=tria->SurfaceAverageVelMisfit();
    3362                 delete tria->material; delete tria;
    3363                 return J;
    3364         }
    3365 }
    3366 /*}}}*/
    3367 IssmDouble Penta::SurfaceLogVxVyMisfit(void){/*{{{*/
    3368 
    3369         IssmDouble J;
    3370         Tria* tria=NULL;
    3371 
    3372         /*inputs: */
    3373         int  approximation;
    3374 
    3375         /*retrieve inputs :*/
    3376         inputs->GetInputValue(&approximation,ApproximationEnum);
    3377 
    3378         /*If on water, return 0: */
    3379         if(!IsIceInElement())return 0;
    3380 
    3381         /*Bail out if this element if:
    3382          * -> Non SSA and not on the surface
    3383          * -> SSA (2d model) and not on bed) */
    3384         if ((approximation!=SSAApproximationEnum && !IsOnSurface()) || (approximation==SSAApproximationEnum && !IsOnBase())){
    3385                 return 0;
    3386         }
    3387         else if (approximation==SSAApproximationEnum){
    3388 
    3389                 /*This element should be collapsed into a tria element at its base. Create this tria element,
    3390                  * and compute SurfaceLogVxVyMisfit*/
    3391                 tria=(Tria*)SpawnTria(0,1,2);
    3392                 J=tria->SurfaceLogVxVyMisfit();
    3393                 delete tria->material; delete tria;
    3394                 return J;
    3395         }
    3396         else{
    3397 
    3398                 tria=(Tria*)SpawnTria(3,4,5);
    3399                 J=tria->SurfaceLogVxVyMisfit();
    3400                 delete tria->material; delete tria;
    3401                 return J;
    3402         }
    3403 }
    3404 /*}}}*/
    3405 IssmDouble Penta::ThicknessAbsGradient(void){/*{{{*/
    3406 
    3407         _error_("Not implemented yet");
    3408 }
    3409 /*}}}*/
    3410 IssmDouble Penta::ThicknessAbsMisfit(void){/*{{{*/
    3411 
    3412         int    approximation;
    3413         IssmDouble J;
    3414         Tria*  tria=NULL;
    3415 
    3416         /*retrieve inputs :*/
    3417         inputs->GetInputValue(&approximation,ApproximationEnum);
    3418 
    3419         /*If on water, return 0: */
    3420         if(!IsIceInElement())return 0;
    3421         _error_("Not implemented yet");
    3422 
    3423         tria=(Tria*)SpawnTria(0,1,2);
    3424         J=tria->ThicknessAbsMisfit();
    3425         delete tria->material; delete tria;
    3426         return J;
    3427 }
    3428 /*}}}*/
    3429 IssmDouble Penta::RheologyBbarAbsGradient(void){/*{{{*/
    3430 
    3431         IssmDouble J;
    3432         Tria*  tria=NULL;
    3433 
    3434         /*If on water, on shelf or not on bed, skip: */
    3435         if(!IsIceInElement() || !IsOnBase()) return 0;
    3436 
    3437         tria=(Tria*)SpawnTria(0,1,2);
    3438         J=tria->RheologyBbarAbsGradient();
    3439         delete tria->material; delete tria;
    3440         return J;
    3441 }
    3442 /*}}}*/
    34433331void       Penta::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data){/*{{{*/
    34443332

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

    r17975 r17976  
    146146                void   ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index);
    147147                void   ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
    148                 IssmDouble RheologyBbarAbsGradient(void);
    149                 IssmDouble ThicknessAbsMisfit(void);
    150                 IssmDouble SurfaceLogVxVyMisfit(void);
    151                 IssmDouble SurfaceAverageVelMisfit(void);
    152                 IssmDouble ThicknessAbsGradient(void);
    153                 IssmDouble ThicknessAlongGradient(void){_error_("not supported");};
    154                 IssmDouble ThicknessAcrossGradient(void){_error_("not supported");};
    155                 IssmDouble BalancethicknessMisfit(void){_error_("not supported");};
    156148                void   InputControlUpdate(IssmDouble scalar,bool save_parameter);
    157149                IssmDouble Misfit(int modelenum,int observationenum,int weightsenum){_error_("not implemented yet");};

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

    r17975 r17976  
    187187                void       ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");};
    188188                void       ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){_error_("not implemented yet");};
    189                 IssmDouble RheologyBbarAbsGradient(void){_error_("not implemented yet");};
    190                 IssmDouble ThicknessAbsMisfit(void){_error_("not implemented yet");};
    191                 IssmDouble ThicknessAbsGradient(void){_error_("not implemented yet");};
    192                 IssmDouble ThicknessAlongGradient(void){_error_("not implemented yet");};
    193                 IssmDouble ThicknessAcrossGradient(void){_error_("not implemented yet");};
    194                 IssmDouble BalancethicknessMisfit(void){_error_("not implemented yet");};
    195                 IssmDouble SurfaceLogVxVyMisfit(void){_error_("not implemented yet");};
    196                 IssmDouble SurfaceAverageVelMisfit(void){_error_("not implemented yet");};
    197189                void       InputControlUpdate(IssmDouble scalar,bool save_parameter){_error_("not implemented yet");};
    198190

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

    r17975 r17976  
    192192                void       ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");};
    193193                void       ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){_error_("not implemented yet");};
    194                 IssmDouble RheologyBbarAbsGradient(void){_error_("not implemented yet");};
    195                 IssmDouble ThicknessAbsMisfit(void){_error_("not implemented yet");};
    196                 IssmDouble ThicknessAbsGradient(void){_error_("not implemented yet");};
    197                 IssmDouble ThicknessAlongGradient(void){_error_("not implemented yet");};
    198                 IssmDouble ThicknessAcrossGradient(void){_error_("not implemented yet");};
    199                 IssmDouble BalancethicknessMisfit(void){_error_("not implemented yet");};
    200                 IssmDouble SurfaceLogVxVyMisfit(void){_error_("not implemented yet");};
    201                 IssmDouble SurfaceAverageVelMisfit(void){_error_("not implemented yet");};
    202194                void       InputControlUpdate(IssmDouble scalar,bool save_parameter){_error_("not implemented yet");};
    203195

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

    r17975 r17976  
    27082708#endif
    27092709
    2710 IssmDouble Tria::BalancethicknessMisfit(void){/*{{{*/
    2711 
    2712         /* Intermediaries */
    2713         IssmDouble Jelem = 0;
    2714         IssmDouble weight;
    2715         IssmDouble Jdet,temp;
    2716         IssmDouble xyz_list[NUMVERTICES][3];
    2717         IssmDouble dH[2];
    2718         IssmDouble vx,vy,H;
    2719         IssmDouble dvx[2],dvy[2];
    2720         IssmDouble dhdt,basal_melting,surface_mass_balance;
    2721         GaussTria *gauss = NULL;
    2722 
    2723         /*If on water, return 0: */
    2724         if(!IsIceInElement()) return 0;
    2725 
    2726         /*Retrieve all inputs we will be needing: */
    2727         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    2728         Input* weights_input              = inputs->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
    2729         Input* thickness_input            = inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    2730         Input* vx_input                   = inputs->GetInput(VxEnum);                                 _assert_(vx_input);
    2731         Input* vy_input                   = inputs->GetInput(VyEnum);                                 _assert_(vy_input);
    2732         Input* surface_mass_balance_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum);         _assert_(surface_mass_balance_input);
    2733         Input* basal_melting_input        = inputs->GetInput(BasalforcingsMeltingRateEnum);           _assert_(basal_melting_input);
    2734         Input* dhdt_input                 = inputs->GetInput(BalancethicknessThickeningRateEnum);     _assert_(dhdt_input);
    2735 
    2736         /* Start looping on the number of gaussian points: */
    2737         gauss=new GaussTria(2);
    2738         for(int ig=gauss->begin();ig<gauss->end();ig++){
    2739 
    2740                 gauss->GaussPoint(ig);
    2741 
    2742                 /* Get Jacobian determinant: */
    2743                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    2744 
    2745                 /*Get all parameters at gaussian point*/
    2746                 weights_input->GetInputValue(&weight,gauss,BalancethicknessMisfitEnum);
    2747                 thickness_input->GetInputValue(&H, gauss);
    2748                 thickness_input->GetInputDerivativeValue(&dH[0],&xyz_list[0][0],gauss);
    2749                 surface_mass_balance_input->GetInputValue(&surface_mass_balance,gauss);
    2750                 basal_melting_input->GetInputValue(&basal_melting,gauss);
    2751                 dhdt_input->GetInputValue(&dhdt,gauss);
    2752                 vx_input->GetInputValue(&vx,gauss);
    2753                 vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
    2754                 vy_input->GetInputValue(&vy,gauss);
    2755                 vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
    2756 
    2757                 /*Balance thickness soft constraint J = 1/2 (div(Hv)-a)^2*/
    2758                 temp  = vx*dH[0]+vy*dH[1]+H*(dvx[0]+dvy[1]) - (surface_mass_balance-basal_melting-dhdt);
    2759                 Jelem+=weight*1/2*temp*temp*Jdet*gauss->weight;
    2760         }
    2761 
    2762         /*Clean up and return*/
    2763         delete gauss;
    2764         return Jelem;
    2765 }
    2766 /*}}}*/
    27672710void       Tria::InputControlUpdate(IssmDouble scalar,bool save_parameter){/*{{{*/
    27682711
    … …  
    34743417}
    34753418/*}}}*/
    3476 IssmDouble Tria::RheologyBbarAbsGradient(void){/*{{{*/
    3477 
    3478         /* Intermediaries */
    3479         IssmDouble Jelem = 0;
    3480         IssmDouble weight;
    3481         IssmDouble Jdet;
    3482         IssmDouble xyz_list[NUMVERTICES][3];
    3483         IssmDouble dp[NDOF2];
    3484         GaussTria *gauss = NULL;
    3485 
    3486         /*retrieve parameters and inputs*/
    3487 
    3488         /*If on water, return 0: */
    3489         if(!IsIceInElement()) return 0;
    3490 
    3491         /*Retrieve all inputs we will be needing: */
    3492         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3493         Input* weights_input  =inputs->GetInput(InversionCostFunctionsCoefficientsEnum);              _assert_(weights_input);
    3494         Input* rheologyb_input=inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
    3495 
    3496         /* Start looping on the number of gaussian points: */
    3497         gauss=new GaussTria(2);
    3498         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3499 
    3500                 gauss->GaussPoint(ig);
    3501 
    3502                 /* Get Jacobian determinant: */
    3503                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3504 
    3505                 /*Get all parameters at gaussian point*/
    3506                 weights_input->GetInputValue(&weight,gauss,RheologyBbarAbsGradientEnum);
    3507                 rheologyb_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
    3508 
    3509                 /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
    3510                 Jelem+=weight*1/2*(dp[0]*dp[0] + dp[1]*dp[1])*Jdet*gauss->weight;
    3511         }
    3512 
    3513         /*Clean up and return*/
    3514         delete gauss;
    3515         return Jelem;
    3516 }
    3517 /*}}}*/
    3518 IssmDouble Tria::SurfaceAverageVelMisfit(void){/*{{{*/
    3519 
    3520         IssmDouble Jelem=0,S,Jdet;
    3521         IssmDouble misfit;
    3522         IssmDouble vx,vy,vxobs,vyobs,weight;
    3523         IssmDouble xyz_list[NUMVERTICES][3];
    3524         GaussTria *gauss=NULL;
    3525 
    3526         /*If on water, return 0: */
    3527         if(!IsIceInElement())return 0;
    3528 
    3529         /* Get node coordinates and dof list: */
    3530         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3531 
    3532         /*Retrieve all inputs we will be needing: */
    3533         inputs->GetInputValue(&S,SurfaceAreaEnum);
    3534         Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
    3535         Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
    3536         Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
    3537         Input* vxobs_input  =inputs->GetInput(InversionVxObsEnum);     _assert_(vxobs_input);
    3538         Input* vyobs_input  =inputs->GetInput(InversionVyObsEnum);     _assert_(vyobs_input);
    3539 
    3540         /* Start  looping on the number of gaussian points: */
    3541         gauss=new GaussTria(3);
    3542         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3543 
    3544                 gauss->GaussPoint(ig);
    3545 
    3546                 /* Get Jacobian determinant: */
    3547                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3548 
    3549                 /*Get all parameters at gaussian point*/
    3550                 weights_input->GetInputValue(&weight,gauss,SurfaceAverageVelMisfitEnum);
    3551                 vx_input->GetInputValue(&vx,gauss);
    3552                 vy_input->GetInputValue(&vy,gauss);
    3553                 vxobs_input->GetInputValue(&vxobs,gauss);
    3554                 vyobs_input->GetInputValue(&vyobs,gauss);
    3555 
    3556                 /*Compute SurfaceAverageVelMisfitEnum:
    3557                  *
    3558                  *      1                    2              2
    3559                  * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
    3560                  *      S                obs            obs
    3561                  */
    3562                 misfit=1/S*sqrt( pow(vx-vxobs,2) + pow(vy-vyobs,2));
    3563 
    3564                 /*Add to cost function*/
    3565                 Jelem+=misfit*weight*Jdet*gauss->weight;
    3566         }
    3567 
    3568         /*clean-up and Return: */
    3569         delete gauss;
    3570         return Jelem;
    3571 }
    3572 /*}}}*/
    3573 IssmDouble Tria::SurfaceLogVxVyMisfit(void){/*{{{*/
    3574 
    3575         IssmDouble Jelem=0, S=0;
    3576         IssmDouble epsvel=2.220446049250313e-16;
    3577         IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/
    3578         IssmDouble misfit, Jdet;
    3579         IssmDouble vx,vy,vxobs,vyobs,weight;
    3580         IssmDouble xyz_list[NUMVERTICES][3];
    3581         GaussTria *gauss=NULL;
    3582 
    3583         /*If on water, return 0: */
    3584         if(!IsIceInElement())return 0;
    3585 
    3586         /* Get node coordinates and dof list: */
    3587         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3588 
    3589         /*Retrieve all inputs we will be needing: */
    3590         Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
    3591         Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
    3592         Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
    3593         Input* vxobs_input  =inputs->GetInput(InversionVxObsEnum);     _assert_(vxobs_input);
    3594         Input* vyobs_input  =inputs->GetInput(InversionVyObsEnum);     _assert_(vyobs_input);
    3595 
    3596         /* Start  looping on the number of gaussian points: */
    3597         gauss=new GaussTria(4);
    3598         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3599 
    3600                 gauss->GaussPoint(ig);
    3601 
    3602                 /* Get Jacobian determinant: */
    3603                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3604 
    3605                 /*Get all parameters at gaussian point*/
    3606                 weights_input->GetInputValue(&weight,gauss,SurfaceLogVxVyMisfitEnum);
    3607                 vx_input->GetInputValue(&vx,gauss);
    3608                 vy_input->GetInputValue(&vy,gauss);
    3609                 vxobs_input->GetInputValue(&vxobs,gauss);
    3610                 vyobs_input->GetInputValue(&vyobs,gauss);
    3611 
    3612                 /*Compute SurfaceRelVelMisfit:
    3613                  *
    3614                  *      1            [        |u| + eps     2          |v| + eps     2  ]
    3615                  * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   | 
    3616                  *      2            [       |u    |+ eps              |v    |+ eps     ]
    3617                  *                              obs                       obs
    3618                  */
    3619                 misfit=0.5*pow(meanvel,2)*(
    3620                                         pow(log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)),2) +
    3621                                         pow(log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)),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 /*}}}*/
    3632 IssmDouble Tria::ThicknessAbsGradient(void){/*{{{*/
    3633 
    3634         /* Intermediaries */
    3635         IssmDouble Jelem = 0;
    3636         IssmDouble weight;
    3637         IssmDouble Jdet;
    3638         IssmDouble xyz_list[NUMVERTICES][3];
    3639         IssmDouble dp[NDOF2];
    3640         GaussTria *gauss = NULL;
    3641 
    3642         /*retrieve parameters and inputs*/
    3643 
    3644         /*If on water, return 0: */
    3645         if(!IsIceInElement()) return 0;
    3646 
    3647         /*Retrieve all inputs we will be needing: */
    3648         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3649         Input* weights_input  =inputs->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
    3650         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    3651 
    3652         /* Start looping on the number of gaussian points: */
    3653         gauss=new GaussTria(2);
    3654         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3655 
    3656                 gauss->GaussPoint(ig);
    3657 
    3658                 /* Get Jacobian determinant: */
    3659                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3660 
    3661                 /*Get all parameters at gaussian point*/
    3662                 weights_input->GetInputValue(&weight,gauss,ThicknessAbsGradientEnum);
    3663                 thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
    3664 
    3665                 /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
    3666                 Jelem+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight;
    3667         }
    3668 
    3669         /*Clean up and return*/
    3670         delete gauss;
    3671         return Jelem;
    3672 }
    3673 /*}}}*/
    3674 IssmDouble Tria::ThicknessAlongGradient(void){/*{{{*/
    3675 
    3676         /* Intermediaries */
    3677         IssmDouble  Jelem = 0;
    3678         IssmDouble  weight;
    3679         IssmDouble  Jdet;
    3680         IssmDouble  xyz_list[NUMVERTICES][3];
    3681         IssmDouble  dp[NDOF2];
    3682         IssmDouble  vx,vy,vel;
    3683         GaussTria  *gauss                    = NULL;
    3684 
    3685         /*retrieve parameters and inputs*/
    3686 
    3687         /*If on water, return 0: */
    3688         if(!IsIceInElement()) return 0;
    3689 
    3690         /*Retrieve all inputs we will be needing: */
    3691         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3692         Input* weights_input  = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    3693         Input* thickness_input= inputs->GetInput(ThicknessEnum);                          _assert_(thickness_input);
    3694         Input* vx_input       = inputs->GetInput(VxEnum);                                 _assert_(vx_input);
    3695         Input* vy_input       = inputs->GetInput(VyEnum);                                 _assert_(vy_input);
    3696 
    3697         /* Start looping on the number of gaussian points: */
    3698         gauss=new GaussTria(2);
    3699         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3700 
    3701                 gauss->GaussPoint(ig);
    3702 
    3703                 /* Get Jacobian determinant: */
    3704                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3705 
    3706                 /*Get all parameters at gaussian point*/
    3707                 weights_input->GetInputValue(&weight,gauss,ThicknessAlongGradientEnum);
    3708                 thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
    3709                 vx_input->GetInputValue(&vx,gauss);
    3710                 vy_input->GetInputValue(&vy,gauss);
    3711                 vel = sqrt(vx*vx+vy*vy);
    3712                 vx  = vx/(vel+1.e-9);
    3713                 vy  = vy/(vel+1.e-9);
    3714 
    3715                 /*J = 1/2 ( vx*dH/dx + vy*dH/dy )^2 */
    3716                 Jelem+=weight*1/2*(vx*dp[0] + vy*dp[1])*(vx*dp[0] + vy*dp[1])*Jdet*gauss->weight;
    3717         }
    3718 
    3719         /*Clean up and return*/
    3720         delete gauss;
    3721         return Jelem;
    3722 }
    3723 /*}}}*/
    3724 IssmDouble Tria::ThicknessAcrossGradient(void){/*{{{*/
    3725 
    3726         /* Intermediaries */
    3727         IssmDouble  Jelem = 0;
    3728         IssmDouble  weight;
    3729         IssmDouble  Jdet;
    3730         IssmDouble  xyz_list[NUMVERTICES][3];
    3731         IssmDouble  dp[NDOF2];
    3732         IssmDouble  vx,vy,vel;
    3733         GaussTria  *gauss                    = NULL;
    3734 
    3735         /*retrieve parameters and inputs*/
    3736 
    3737         /*If on water, return 0: */
    3738         if(!IsIceInElement()) return 0;
    3739 
    3740         /*Retrieve all inputs we will be needing: */
    3741         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3742         Input* weights_input  = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    3743         Input* thickness_input= inputs->GetInput(ThicknessEnum);                          _assert_(thickness_input);
    3744         Input* vx_input       = inputs->GetInput(VxEnum);                                 _assert_(vx_input);
    3745         Input* vy_input       = inputs->GetInput(VyEnum);                                 _assert_(vy_input);
    3746 
    3747         /* Start looping on the number of gaussian points: */
    3748         gauss=new GaussTria(2);
    3749         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3750 
    3751                 gauss->GaussPoint(ig);
    3752 
    3753                 /* Get Jacobian determinant: */
    3754                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3755 
    3756                 /*Get all parameters at gaussian point*/
    3757                 weights_input->GetInputValue(&weight,gauss,ThicknessAcrossGradientEnum);
    3758                 thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
    3759                 vx_input->GetInputValue(&vx,gauss);
    3760                 vy_input->GetInputValue(&vy,gauss);
    3761                 vel = sqrt(vx*vx+vy*vy);
    3762                 vx  = vx/(vel+1.e-9);
    3763                 vy  = vy/(vel+1.e-9);
    3764 
    3765                 /*J = 1/2 ( -vy*dH/dx + vx*dH/dy )^2 */
    3766                 Jelem+=weight*1/2*(-vy*dp[0] + vx*dp[1])*(-vy*dp[0] + vx*dp[1])*Jdet*gauss->weight;
    3767         }
    3768 
    3769         /*Clean up and return*/
    3770         delete gauss;
    3771         return Jelem;
    3772 }
    3773 /*}}}*/
    3774 IssmDouble Tria::ThicknessAbsMisfit(void){/*{{{*/
    3775 
    3776         /*Intermediaries*/
    3777         IssmDouble thickness,thicknessobs,weight;
    3778         IssmDouble Jdet;
    3779         IssmDouble Jelem = 0;
    3780         IssmDouble xyz_list[NUMVERTICES][3];
    3781         GaussTria *gauss = NULL;
    3782         IssmDouble dH[2];
    3783 
    3784         /*If on water, return 0: */
    3785         if(!IsIceInElement())return 0;
    3786 
    3787         /*Retrieve all inputs we will be needing: */
    3788         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3789         Input* thickness_input   =inputs->GetInput(ThicknessEnum);   _assert_(thickness_input);
    3790         Input* thicknessobs_input=inputs->GetInput(InversionThicknessObsEnum);_assert_(thicknessobs_input);
    3791         Input* weights_input     =inputs->GetInput(InversionCostFunctionsCoefficientsEnum);     _assert_(weights_input);
    3792 
    3793         /* Start  looping on the number of gaussian points: */
    3794         gauss=new GaussTria(2);
    3795         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3796 
    3797                 gauss->GaussPoint(ig);
    3798 
    3799                 /* Get Jacobian determinant: */
    3800                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3801 
    3802                 /*Get parameters at gauss point*/
    3803                 thickness_input->GetInputValue(&thickness,gauss);
    3804                 thickness_input->GetInputDerivativeValue(&dH[0],&xyz_list[0][0],gauss);
    3805                 thicknessobs_input->GetInputValue(&thicknessobs,gauss);
    3806                 weights_input->GetInputValue(&weight,gauss,ThicknessAbsMisfitEnum);
    3807 
    3808                 /*compute ThicknessAbsMisfit*/
    3809                 Jelem+=0.5*(thickness-thicknessobs)*(thickness-thicknessobs)*weight*Jdet*gauss->weight;
    3810         }
    3811 
    3812         /* clean up and Return: */
    3813         delete gauss;
    3814         return Jelem;
    3815 }
    3816 /*}}}*/
    38173419void       Tria::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data){/*{{{*/
    38183420

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

    r17975 r17976  
    156156                void       ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index);
    157157                void       ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
    158                 IssmDouble RheologyBbarAbsGradient(void);
    159                 IssmDouble ThicknessAbsMisfit(void);
    160                 IssmDouble ThicknessAbsGradient(void);
    161                 IssmDouble ThicknessAlongGradient(void);
    162                 IssmDouble ThicknessAcrossGradient(void);
    163                 IssmDouble BalancethicknessMisfit(void);
    164                 IssmDouble SurfaceLogVxVyMisfit(void);
    165                 IssmDouble SurfaceAverageVelMisfit(void);
    166158                void       InputControlUpdate(IssmDouble scalar,bool save_parameter);
    167159

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r17907 r17976  
    11341134void FemModel::BalancethicknessMisfitx(IssmDouble* presponse){/*{{{*/
    11351135
    1136         IssmDouble J = 0;
     1136        /*output: */
     1137        IssmDouble J=0.;
    11371138        IssmDouble J_sum;
    11381139
    1139         for(int i=0;i<this->elements->Size();i++){
    1140                 Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
    1141                 J+=element->BalancethicknessMisfit();
    1142         }
     1140        IssmDouble  weight,vx,vy,H,dvx[2],dvy[2],dH[2];
     1141        IssmDouble  temp,Jdet,dhdt,basal_melting,surface_mass_balance;
     1142        IssmDouble* xyz_list = NULL;
     1143        IssmDouble  dp[3];
     1144
     1145        /*Compute Misfit: */
     1146        for(int i=0;i<elements->Size();i++){
     1147                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
     1148
     1149                /*If on water, return 0: */
     1150                if(!element->IsIceInElement()) continue;
     1151
     1152                /* Get node coordinates*/
     1153                element->GetVerticesCoordinates(&xyz_list);
     1154                Input* weights_input              = element->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
     1155                Input* thickness_input            = element->GetInput(ThicknessEnum); _assert_(thickness_input);
     1156                Input* vx_input                   = element->GetInput(VxEnum);                                 _assert_(vx_input);
     1157                Input* vy_input                   = element->GetInput(VyEnum);                                 _assert_(vy_input);
     1158                Input* surface_mass_balance_input = element->GetInput(SurfaceforcingsMassBalanceEnum);         _assert_(surface_mass_balance_input);
     1159                Input* basal_melting_input        = element->GetInput(BasalforcingsMeltingRateEnum);           _assert_(basal_melting_input);
     1160                Input* dhdt_input                 = element->GetInput(BalancethicknessThickeningRateEnum);     _assert_(dhdt_input);
     1161
     1162                /* Start  looping on the number of gaussian points: */
     1163                Gauss* gauss=element->NewGauss(2);
     1164                for(int ig=gauss->begin();ig<gauss->end();ig++){
     1165
     1166                        gauss->GaussPoint(ig);
     1167
     1168                        /* Get Jacobian determinant: */
     1169                        element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     1170
     1171                        /*Get all parameters at gaussian point*/
     1172                        weights_input->GetInputValue(&weight,gauss,BalancethicknessMisfitEnum);
     1173                        thickness_input->GetInputValue(&H, gauss);
     1174                        thickness_input->GetInputDerivativeValue(&dH[0],xyz_list,gauss);
     1175                        surface_mass_balance_input->GetInputValue(&surface_mass_balance,gauss);
     1176                        basal_melting_input->GetInputValue(&basal_melting,gauss);
     1177                        dhdt_input->GetInputValue(&dhdt,gauss);
     1178                        vx_input->GetInputValue(&vx,gauss);
     1179                        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
     1180                        vy_input->GetInputValue(&vy,gauss);
     1181                        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
     1182
     1183                        /*Balance thickness soft constraint J = 1/2 (div(Hv)-a)^2*/
     1184                        temp  = vx*dH[0]+vy*dH[1]+H*(dvx[0]+dvy[1]) - (surface_mass_balance-basal_melting-dhdt);
     1185                        J    +=weight*1/2*temp*temp*Jdet*gauss->weight;
     1186                }
     1187
     1188                /*clean up and Return: */
     1189                xDelete<IssmDouble>(xyz_list);
     1190                delete gauss;
     1191        }
     1192
     1193        /*Sum all J from all cpus of the cluster:*/
    11431194        ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    11441195        ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    … …  
    11511202void FemModel::ThicknessAbsGradientx( IssmDouble* pJ){/*{{{*/
    11521203
    1153         /*Intermediary*/
    1154         int i;
    1155         Element* element=NULL;
    1156 
    11571204        /*output: */
    1158         IssmDouble J=0;
     1205        IssmDouble J=0.;
    11591206        IssmDouble J_sum;
    11601207
     1208        IssmDouble  thickness,thicknessobs,weight;
     1209        IssmDouble  Jdet;
     1210        IssmDouble* xyz_list = NULL;
     1211        IssmDouble  dp[3];
     1212
    11611213        /*Compute Misfit: */
    1162         for (i=0;i<elements->Size();i++){
    1163                 element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
    1164                 J+=element->ThicknessAbsGradient();
     1214        for(int i=0;i<elements->Size();i++){
     1215                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
     1216
     1217                /*If on water, return 0: */
     1218                if(!element->IsIceInElement()) continue;
     1219
     1220                /* Get node coordinates*/
     1221                element->GetVerticesCoordinates(&xyz_list);
     1222
     1223                /*Retrieve all inputs we will be needing: */
     1224                Input* weights_input   =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     1225                Input* thickness_input =element->GetInput(ThicknessEnum);                          _assert_(thickness_input);
     1226
     1227                /* Start  looping on the number of gaussian points: */
     1228                Gauss* gauss=element->NewGauss(2);
     1229                for(int ig=gauss->begin();ig<gauss->end();ig++){
     1230
     1231                        gauss->GaussPoint(ig);
     1232
     1233                        /* Get Jacobian determinant: */
     1234                        element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     1235
     1236                        /*Get all parameters at gaussian point*/
     1237                        weights_input->GetInputValue(&weight,gauss,ThicknessAcrossGradientEnum);
     1238                        thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
     1239
     1240                        /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
     1241                        J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight;
     1242                }
     1243
     1244                /*clean up and Return: */
     1245                xDelete<IssmDouble>(xyz_list);
     1246                delete gauss;
    11651247        }
    11661248

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

    r17972 r17976  
    3232
    3333        int         domaintype,numcomponents;
    34         bool        surfaceintegral;
    3534        IssmDouble  Jelem=0.;
    3635        IssmDouble  misfit,Jdet;
    … …  
    4746        element->FindParam(&domaintype,DomainTypeEnum);
    4847        switch(domaintype){
    49                 case Domain2DverticalEnum:
    50                         surfaceintegral = true;
    51                         numcomponents   = 1;
    52                         break;
    53                 case Domain3DEnum:
    54                         surfaceintegral = true;
    55                         numcomponents   = 2;
    56                         break;
    57                 case Domain2DhorizontalEnum:
    58                         surfaceintegral = false;
    59                         numcomponents   = 2;
    60                         break;
     48                case Domain2DverticalEnum:   numcomponents   = 1; break;
     49                case Domain3DEnum:           numcomponents   = 2; break;
     50                case Domain2DhorizontalEnum: numcomponents   = 2; break;
    6151                default: _error_("not supported yet");
    6252        }

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

    r16314 r17976  
    1010void RheologyBbarAbsGradientx( 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->RheologyBbarAbsGradient();
     17        for(int i=0;i<elements->Size();i++){
     18                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
     19                J+=RheologyBbarAbsGradient(element);
    2420        }
    2521
    … …  
    3228        *pJ=J;
    3329}
     30
     31IssmDouble RheologyBbarAbsGradient(Element* element){
     32
     33        int         domaintype,numcomponents;
     34        IssmDouble  Jelem=0.;
     35        IssmDouble  misfit,Jdet;
     36        IssmDouble  dp[2],weight;
     37        IssmDouble* xyz_list      = NULL;
     38
     39        /*Get basal element*/
     40        if(!element->IsOnBase()) return 0.;
     41
     42        /*If on water, return 0: */
     43        if(!element->IsIceInElement()) return 0.;
     44
     45        /*Get problem dimension*/
     46        element->FindParam(&domaintype,DomainTypeEnum);
     47        switch(domaintype){
     48                case Domain2DverticalEnum:   numcomponents   = 1; break;
     49                case Domain3DEnum:           numcomponents   = 2; break;
     50                case Domain2DhorizontalEnum: numcomponents   = 2; break;
     51                default: _error_("not supported yet");
     52        }
     53
     54        /*Spawn basal element*/
     55        Element* basalelement = element->SpawnBasalElement();
     56
     57        /* Get node coordinates*/
     58        basalelement->GetVerticesCoordinates(&xyz_list);
     59
     60        /*Retrieve all inputs we will be needing: */
     61        Input* weights_input=basalelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     62        Input* rheologyb_input=basalelement->GetInput(MaterialsRheologyBbarEnum);                  _assert_(rheologyb_input);
     63
     64        /* Start  looping on the number of gaussian points: */
     65        Gauss* gauss=basalelement->NewGauss(2);
     66        for(int ig=gauss->begin();ig<gauss->end();ig++){
     67
     68                gauss->GaussPoint(ig);
     69
     70                /* Get Jacobian determinant: */
     71                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     72
     73                /*Get all parameters at gaussian point*/
     74                weights_input->GetInputValue(&weight,gauss,RheologyBbarAbsGradientEnum);
     75                rheologyb_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
     76
     77                /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
     78                Jelem+=weight*1/2*(dp[0]*dp[0] + dp[1]*dp[1])*Jdet*gauss->weight;
     79        }
     80
     81        /*clean up and Return: */
     82        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     83        xDelete<IssmDouble>(xyz_list);
     84        delete gauss;
     85        return Jelem;
     86}

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

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

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

    r17972 r17976  
    3232
    3333        int        domaintype,numcomponents;
    34         bool       surfaceintegral;
    3534        IssmDouble Jelem=0.;
    3635        IssmDouble misfit,Jdet;
    … …  
    4746        element->FindParam(&domaintype,DomainTypeEnum);
    4847        switch(domaintype){
    49                 case Domain2DverticalEnum:
    50                         surfaceintegral = true;
    51                         numcomponents   = 1;
    52                         break;
    53                 case Domain3DEnum:
    54                         surfaceintegral = true;
    55                         numcomponents   = 2;
    56                         break;
    57                 case Domain2DhorizontalEnum:
    58                         surfaceintegral = false;
    59                         numcomponents   = 2;
    60                         break;
     48                case Domain2DverticalEnum:   numcomponents   = 1; break;
     49                case Domain3DEnum:           numcomponents   = 2; break;
     50                case Domain2DhorizontalEnum: numcomponents   = 2; break;
    6151                default: _error_("not supported yet");
    6252        }

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

    r16314 r17976  
    1111void SurfaceAverageVelMisfitx(IssmDouble* pJ,FemModel* femmodel){
    1212
    13         /*Intermediary*/
    14         Element* element=NULL;
    15 
    1613        /*output: */
    17         IssmDouble J = 0.;
     14        IssmDouble J=0.;
    1815        IssmDouble J_sum;
    19 
    20         /*Compute surface area and add to elements inputs */
    21         SurfaceAreax(NULL,femmodel);
    2216
    2317        /*Compute Misfit: */
    2418        for(int i=0;i<femmodel->elements->Size();i++){
    25                 element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
    26                 J+=element->SurfaceAverageVelMisfit();
     19                Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
     20                J+=SurfaceAverageVelMisfit(element);
    2721        }
    2822
    … …  
    3529        *pJ=J;
    3630}
     31
     32IssmDouble SurfaceAverageVelMisfit(Element* element){
     33
     34        int        domaintype,numcomponents;
     35        IssmDouble Jelem=0.;
     36        IssmDouble misfit,S,Jdet;
     37        IssmDouble vx,vy,vxobs,vyobs,weight;
     38        IssmDouble* xyz_list = NULL;
     39
     40        /*Get basal element*/
     41        if(!element->IsOnSurface()) return 0.;
     42
     43        /*If on water, return 0: */
     44        if(!element->IsIceInElement()) return 0.;
     45
     46        /*Get problem dimension*/
     47        element->FindParam(&domaintype,DomainTypeEnum);
     48        switch(domaintype){
     49                case Domain2DverticalEnum:
     50                        numcomponents   = 1;
     51                        break;
     52                case Domain3DEnum:
     53                        numcomponents   = 2;
     54                        break;
     55                case Domain2DhorizontalEnum:
     56                        numcomponents   = 2;
     57                        break;
     58                default: _error_("not supported yet");
     59        }
     60
     61        /*Spawn surface element*/
     62        Element* topelement = element->SpawnTopElement();
     63
     64        /* Get node coordinates*/
     65        topelement->GetVerticesCoordinates(&xyz_list);
     66
     67        /*Retrieve all inputs we will be needing: */
     68        topelement->GetInputValue(&S,SurfaceAreaEnum);
     69        Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     70        Input* vx_input     =topelement->GetInput(VxEnum);                                 _assert_(vx_input);
     71        Input* vxobs_input  =topelement->GetInput(InversionVxObsEnum);                     _assert_(vxobs_input);
     72        Input* vy_input     = NULL;
     73        Input* vyobs_input  = NULL;
     74        if(numcomponents==2){
     75                vy_input    =topelement->GetInput(VyEnum);              _assert_(vy_input);
     76                vyobs_input =topelement->GetInput(InversionVyObsEnum);  _assert_(vyobs_input);
     77        }
     78
     79        /* Start  looping on the number of gaussian points: */
     80        Gauss* gauss=topelement->NewGauss(3);
     81        for(int ig=gauss->begin();ig<gauss->end();ig++){
     82
     83                gauss->GaussPoint(ig);
     84
     85                /* Get Jacobian determinant: */
     86                topelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     87
     88                /*Get all parameters at gaussian point*/
     89                weights_input->GetInputValue(&weight,gauss,SurfaceAverageVelMisfitEnum);
     90                vx_input->GetInputValue(&vx,gauss);
     91                vxobs_input->GetInputValue(&vxobs,gauss);
     92                if(numcomponents==2){
     93                        vy_input->GetInputValue(&vy,gauss);
     94                        vyobs_input->GetInputValue(&vyobs,gauss);
     95                }
     96
     97                /*Compute SurfaceAverageVelMisfitEnum:
     98                 *
     99                 *      1                    2              2
     100                 * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
     101                 *      S                obs            obs
     102                 */
     103                if(numcomponents==1){
     104                        misfit=1/S*(vx-vxobs)*(vx-vxobs);
     105                }
     106                else{
     107                        misfit=1/S*sqrt( pow(vx-vxobs,2) + pow(vy-vyobs,2));
     108                }
     109
     110                /*Add to cost function*/
     111                Jelem+=misfit*weight*Jdet*gauss->weight;
     112        }
     113
     114        /*clean up and Return: */
     115        if(domaintype!=Domain2DhorizontalEnum){topelement->DeleteMaterials(); delete topelement;};
     116        xDelete<IssmDouble>(xyz_list);
     117        delete gauss;
     118        return Jelem;
     119}

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

    r16314 r17976  
    1010/* local prototypes: */
    1111void SurfaceAverageVelMisfitx(IssmDouble* pJ,FemModel* femmodel);
     12IssmDouble SurfaceAverageVelMisfit(Element* element);
    1213
    1314#endif

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

    r17975 r17976  
    3232
    3333        int        domaintype,numcomponents;
    34         bool       surfaceintegral;
    3534        IssmDouble Jelem=0.;
    3635        IssmDouble epsvel=2.220446049250313e-16;
    … …  
    5049        element->FindParam(&domaintype,DomainTypeEnum);
    5150        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;
     51                case Domain2DverticalEnum:   numcomponents   = 1; break;
     52                case Domain3DEnum:           numcomponents   = 2; break;
     53                case Domain2DhorizontalEnum: numcomponents   = 2; break;
    6454                default: _error_("not supported yet");
    6555        }

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

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

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

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

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

    r17975 r17976  
    3232
    3333        int        domaintype,numcomponents;
    34         bool       surfaceintegral;
    3534        IssmDouble Jelem=0.;
    3635        IssmDouble misfit,Jdet,scalex,scaley;
    … …  
    4948        element->FindParam(&domaintype,DomainTypeEnum);
    5049        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;
     50                case Domain2DverticalEnum:   numcomponents   = 1; break;
     51                case Domain3DEnum:           numcomponents   = 2; break;
     52                case Domain2DhorizontalEnum: numcomponents   = 2; break;
    6353                default: _error_("not supported yet");
    6454        }

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

    r16314 r17976  
    1010void ThicknessAbsMisfitx( 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->ThicknessAbsMisfit();
     17        for(int i=0;i<elements->Size();i++){
     18                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
     19                J+=ThicknessAbsMisfit(element);
    2420        }
    2521
    … …  
    3228        *pJ=J;
    3329}
     30
     31IssmDouble ThicknessAbsMisfit(Element* element){
     32
     33        IssmDouble  thickness,thicknessobs,weight;
     34        IssmDouble  Jelem=0.;
     35        IssmDouble  misfit,Jdet;
     36        IssmDouble* xyz_list = NULL;
     37
     38        /*If on water, return 0: */
     39        if(!element->IsIceInElement()) return 0.;
     40
     41        /* Get node coordinates*/
     42        element->GetVerticesCoordinates(&xyz_list);
     43
     44        /*Retrieve all inputs we will be needing: */
     45        Input* weights_input     =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     46        Input* thickness_input   =element->GetInput(ThicknessEnum);                          _assert_(thickness_input);
     47        Input* thicknessobs_input=element->GetInput(InversionThicknessObsEnum);              _assert_(thicknessobs_input);
     48
     49        /* Start  looping on the number of gaussian points: */
     50        Gauss* gauss=element->NewGauss(2);
     51        for(int ig=gauss->begin();ig<gauss->end();ig++){
     52
     53                gauss->GaussPoint(ig);
     54
     55                /* Get Jacobian determinant: */
     56                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     57
     58                /*Get all parameters at gaussian point*/
     59                weights_input->GetInputValue(&weight,gauss,ThicknessAbsMisfitEnum);
     60                thickness_input->GetInputValue(&thickness,gauss);
     61                thicknessobs_input->GetInputValue(&thicknessobs,gauss);
     62
     63                /*Compute ThicknessAbsMisfitEnum*/
     64                Jelem+=0.5*(thickness-thicknessobs)*(thickness-thicknessobs)*weight*Jdet*gauss->weight;
     65        }
     66
     67        /*clean up and Return: */
     68        xDelete<IssmDouble>(xyz_list);
     69        delete gauss;
     70        return Jelem;
     71}

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

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

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

    r16314 r17976  
    1010void ThicknessAcrossGradientx( 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->ThicknessAcrossGradient();
     17        for(int i=0;i<elements->Size();i++){
     18                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
     19                J+=ThicknessAcrossGradient(element);
    2420        }
    2521
    … …  
    3228        *pJ=J;
    3329}
     30
     31IssmDouble ThicknessAcrossGradient(Element* element){
     32
     33        IssmDouble  thickness,thicknessobs,weight;
     34        IssmDouble  Jelem=0.;
     35        IssmDouble  misfit,Jdet;
     36        IssmDouble* xyz_list = NULL;
     37        IssmDouble  dp[3];
     38        IssmDouble  vx,vy,vel;
     39
     40        /*If on water, return 0: */
     41        if(!element->IsIceInElement()) return 0.;
     42
     43        /* Get node coordinates*/
     44        element->GetVerticesCoordinates(&xyz_list);
     45
     46        /*Retrieve all inputs we will be needing: */
     47        Input* weights_input   =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     48        Input* thickness_input =element->GetInput(ThicknessEnum);                          _assert_(thickness_input);
     49        Input* vx_input        =element->GetInput(VxEnum);                                 _assert_(vx_input);
     50        Input* vy_input        =element->GetInput(VyEnum);                                 _assert_(vy_input);
     51
     52        /* Start  looping on the number of gaussian points: */
     53        Gauss* gauss=element->NewGauss(2);
     54        for(int ig=gauss->begin();ig<gauss->end();ig++){
     55
     56                gauss->GaussPoint(ig);
     57
     58                /* Get Jacobian determinant: */
     59                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     60
     61                /*Get all parameters at gaussian point*/
     62                weights_input->GetInputValue(&weight,gauss,ThicknessAcrossGradientEnum);
     63                thickness_input->GetInputValue(&thickness,gauss);
     64                thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
     65                vx_input->GetInputValue(&vx,gauss);
     66                vy_input->GetInputValue(&vy,gauss);
     67                vel = sqrt(vx*vx+vy*vy);
     68                vx  = vx/(vel+1.e-9);
     69                vy  = vy/(vel+1.e-9);
     70
     71                /*J = 1/2 ( -vy*dH/dx + vx*dH/dy )^2 */
     72                Jelem+=weight*1/2*(-vy*dp[0] + vx*dp[1])*(-vy*dp[0] + vx*dp[1])*Jdet*gauss->weight;
     73        }
     74
     75        /*clean up and Return: */
     76        xDelete<IssmDouble>(xyz_list);
     77        delete gauss;
     78        return Jelem;
     79}

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

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

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

    r16314 r17976  
    1010void ThicknessAlongGradientx( 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->ThicknessAlongGradient();
     17        for(int i=0;i<elements->Size();i++){
     18                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
     19                J+=ThicknessAlongGradient(element);
    2420        }
    2521
    … …  
    3228        *pJ=J;
    3329}
     30
     31IssmDouble ThicknessAlongGradient(Element* element){
     32
     33        IssmDouble  thickness,thicknessobs,weight;
     34        IssmDouble  Jelem=0.;
     35        IssmDouble  misfit,Jdet;
     36        IssmDouble* xyz_list = NULL;
     37        IssmDouble  dp[3];
     38        IssmDouble  vx,vy,vel;
     39
     40        /*If on water, return 0: */
     41        if(!element->IsIceInElement()) return 0.;
     42
     43        /* Get node coordinates*/
     44        element->GetVerticesCoordinates(&xyz_list);
     45
     46        /*Retrieve all inputs we will be needing: */
     47        Input* weights_input   =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     48        Input* thickness_input =element->GetInput(ThicknessEnum);                          _assert_(thickness_input);
     49        Input* vx_input        =element->GetInput(VxEnum);                                 _assert_(vx_input);
     50        Input* vy_input        =element->GetInput(VyEnum);                                 _assert_(vy_input);
     51
     52        /* Start  looping on the number of gaussian points: */
     53        Gauss* gauss=element->NewGauss(2);
     54        for(int ig=gauss->begin();ig<gauss->end();ig++){
     55
     56                gauss->GaussPoint(ig);
     57
     58                /* Get Jacobian determinant: */
     59                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     60
     61                /*Get all parameters at gaussian point*/
     62                weights_input->GetInputValue(&weight,gauss,ThicknessAlongGradientEnum);
     63                thickness_input->GetInputValue(&thickness,gauss);
     64                thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
     65                vx_input->GetInputValue(&vx,gauss);
     66                vy_input->GetInputValue(&vy,gauss);
     67                vel = sqrt(vx*vx+vy*vy);
     68                vx  = vx/(vel+1.e-9);
     69                vy  = vy/(vel+1.e-9);
     70
     71                /*J = 1/2 ( vx*dH/dx + vy*dH/dy )^2 */
     72                Jelem+=weight*1/2*(vx*dp[0] + vy*dp[1])*(vx*dp[0] + vy*dp[1])*Jdet*gauss->weight;
     73        }
     74
     75        /*clean up and Return: */
     76        xDelete<IssmDouble>(xyz_list);
     77        delete gauss;
     78        return Jelem;
     79}

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

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