Changeset 17975
- Timestamp:
- 05/09/14 14:51:29 (11 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 11 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/Elements/Element.h
r17926 r17975 277 277 virtual void Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index)=0; 278 278 virtual IssmDouble ThicknessAbsMisfit(void)=0; 279 virtual IssmDouble SurfaceAbsVelMisfit(void)=0;280 virtual IssmDouble SurfaceRelVelMisfit(void)=0;281 virtual IssmDouble SurfaceLogVelMisfit(void)=0;282 279 virtual IssmDouble SurfaceLogVxVyMisfit(void)=0; 283 280 virtual IssmDouble SurfaceAverageVelMisfit(void)=0; … … 287 284 virtual IssmDouble BalancethicknessMisfit(void)=0; 288 285 virtual IssmDouble RheologyBbarAbsGradient(void)=0; 289 virtual IssmDouble DragCoefficientAbsGradient(void)=0;290 286 virtual void ControlInputGetGradient(Vector<IssmDouble>* gradient,int enum_type,int control_index)=0; 291 287 virtual void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index)=0; -
issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
r17967 r17975 3365 3365 } 3366 3366 /*}}}*/ 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 surface3381 * -> 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 surface3417 * -> 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 /*}}}*/3439 3367 IssmDouble Penta::SurfaceLogVxVyMisfit(void){/*{{{*/ 3440 3368 … … 3475 3403 } 3476 3404 /*}}}*/ 3477 IssmDouble Penta::SurfaceRelVelMisfit(void){/*{{{*/ 3405 IssmDouble Penta::ThicknessAbsGradient(void){/*{{{*/ 3406 3407 _error_("Not implemented yet"); 3408 } 3409 /*}}}*/ 3410 IssmDouble Penta::ThicknessAbsMisfit(void){/*{{{*/ 3478 3411 3479 3412 int approximation; … … 3486 3419 /*If on water, return 0: */ 3487 3420 if(!IsIceInElement())return 0; 3488 3489 /*Bail out if this element if:3490 * -> Non SSA and not on the surface3491 * -> 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;3529 3421 _error_("Not implemented yet"); 3530 3422 3531 3423 tria=(Tria*)SpawnTria(0,1,2); 3532 3424 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 13546 J=tria->DragCoefficientAbsGradient();3547 3425 delete tria->material; delete tria; 3548 3426 return J; -
issm/trunk-jpl/src/c/classes/Elements/Penta.h
r17926 r17975 132 132 #endif 133 133 134 IssmDouble DragCoefficientAbsGradient(void);135 134 void GradientIndexing(int* indexing,int control_index); 136 135 void Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index); … … 149 148 IssmDouble RheologyBbarAbsGradient(void); 150 149 IssmDouble ThicknessAbsMisfit(void); 151 IssmDouble SurfaceAbsVelMisfit(void);152 IssmDouble SurfaceRelVelMisfit(void);153 IssmDouble SurfaceLogVelMisfit(void);154 150 IssmDouble SurfaceLogVxVyMisfit(void); 155 151 IssmDouble SurfaceAverageVelMisfit(void); -
issm/trunk-jpl/src/c/classes/Elements/Seg.h
r17972 r17975 168 168 #endif 169 169 170 IssmDouble DragCoefficientAbsGradient(void){_error_("not implemented yet");};171 170 void GradientIndexing(int* indexing,int control_index); 172 171 void Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index){_error_("not implemented yet");}; … … 190 189 IssmDouble RheologyBbarAbsGradient(void){_error_("not implemented yet");}; 191 190 IssmDouble ThicknessAbsMisfit(void){_error_("not implemented yet");}; 192 IssmDouble SurfaceAbsVelMisfit(void){_error_("not implemented yet");};193 191 IssmDouble ThicknessAbsGradient(void){_error_("not implemented yet");}; 194 192 IssmDouble ThicknessAlongGradient(void){_error_("not implemented yet");}; 195 193 IssmDouble ThicknessAcrossGradient(void){_error_("not implemented yet");}; 196 194 IssmDouble BalancethicknessMisfit(void){_error_("not implemented yet");}; 197 IssmDouble SurfaceRelVelMisfit(void){_error_("not implemented yet");};198 IssmDouble SurfaceLogVelMisfit(void){_error_("not implemented yet");};199 195 IssmDouble SurfaceLogVxVyMisfit(void){_error_("not implemented yet");}; 200 196 IssmDouble SurfaceAverageVelMisfit(void){_error_("not implemented yet");}; -
issm/trunk-jpl/src/c/classes/Elements/Tetra.h
r17926 r17975 194 194 IssmDouble RheologyBbarAbsGradient(void){_error_("not implemented yet");}; 195 195 IssmDouble ThicknessAbsMisfit(void){_error_("not implemented yet");}; 196 IssmDouble SurfaceAbsVelMisfit(void){_error_("not implemented yet");};197 196 IssmDouble ThicknessAbsGradient(void){_error_("not implemented yet");}; 198 197 IssmDouble ThicknessAlongGradient(void){_error_("not implemented yet");}; 199 198 IssmDouble ThicknessAcrossGradient(void){_error_("not implemented yet");}; 200 199 IssmDouble BalancethicknessMisfit(void){_error_("not implemented yet");}; 201 IssmDouble SurfaceRelVelMisfit(void){_error_("not implemented yet");};202 IssmDouble SurfaceLogVelMisfit(void){_error_("not implemented yet");};203 200 IssmDouble SurfaceLogVxVyMisfit(void){_error_("not implemented yet");}; 204 201 IssmDouble SurfaceAverageVelMisfit(void){_error_("not implemented yet");}; -
issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
r17972 r17975 3571 3571 } 3572 3572 /*}}}*/ 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 ] 23615 * J = 4 \bar{v}^2 | log ( ----------- ) |3616 * [ vel + eps ]3617 * obs3618 */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 /*}}}*/3632 3573 IssmDouble Tria::SurfaceLogVxVyMisfit(void){/*{{{*/ 3633 3574 … … 3689 3630 } 3690 3631 /*}}}*/ 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 obs3792 */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 /*}}}*/3806 3632 IssmDouble Tria::ThicknessAbsGradient(void){/*{{{*/ 3807 3633 … … 3985 3811 3986 3812 /* 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*/4029 3813 delete gauss; 4030 3814 return Jelem; -
issm/trunk-jpl/src/c/classes/Elements/Tria.h
r17926 r17975 137 137 #endif 138 138 139 IssmDouble DragCoefficientAbsGradient(void);140 139 void GradientIndexing(int* indexing,int control_index); 141 140 void Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index); … … 159 158 IssmDouble RheologyBbarAbsGradient(void); 160 159 IssmDouble ThicknessAbsMisfit(void); 161 IssmDouble SurfaceAbsVelMisfit(void);162 160 IssmDouble ThicknessAbsGradient(void); 163 161 IssmDouble ThicknessAlongGradient(void); 164 162 IssmDouble ThicknessAcrossGradient(void); 165 163 IssmDouble BalancethicknessMisfit(void); 166 IssmDouble SurfaceRelVelMisfit(void);167 IssmDouble SurfaceLogVelMisfit(void);168 164 IssmDouble SurfaceLogVxVyMisfit(void); 169 165 IssmDouble SurfaceAverageVelMisfit(void); -
issm/trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp
r16314 r17975 10 10 void SurfaceLogVelMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters){ 11 11 12 /*Intermediary*/13 int i;14 Element* element=NULL;15 16 12 /*output: */ 17 IssmDouble J=0 ;13 IssmDouble J=0.; 18 14 IssmDouble J_sum; 19 15 20 16 /*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); 24 20 } 25 21 … … 32 28 *pJ=J; 33 29 } 30 31 IssmDouble 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 10 10 /* local prototypes: */ 11 11 void SurfaceLogVelMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters); 12 IssmDouble SurfaceLogVelMisfit(Element* element); 12 13 13 14 #endif -
issm/trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp
r16314 r17975 10 10 void SurfaceRelVelMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters){ 11 11 12 /*Intermediary*/13 int i;14 Element* element=NULL;15 16 12 /*output: */ 17 IssmDouble J=0 ;13 IssmDouble J=0.; 18 14 IssmDouble J_sum; 19 15 20 16 /*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); 24 20 } 25 21 … … 32 28 *pJ=J; 33 29 } 30 31 IssmDouble 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 10 10 /* local prototypes: */ 11 11 void SurfaceRelVelMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters); 12 IssmDouble SurfaceRelVelMisfit(Element* element); 12 13 13 14 #endif
Note:
See TracChangeset
for help on using the changeset viewer.