Changeset 18062


Ignore:
Timestamp:
05/26/14 20:12:51 (11 years ago)
Author:
Mathieu Morlighem
Message:

DEL: removed all gradients of elements, now done by adjoint

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

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp

    r18057 r18062  
    480480}/*}}}*/
    481481void BalancethicknessAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
    482         _error_("Not implemented yet");
     482        /* WARNING: this gradient is valid for Soft balance thickness only */
     483
     484        /*If on water, grad = 0: */
     485        if(!element->IsIceInElement()) return;
     486
     487        /*Intermediaries*/
     488        IssmDouble Jdet,weight;
     489        IssmDouble thickness,thicknessobs,dH[3],dp[3];
     490        IssmDouble  vx,vy,vel,dvx[2],dvy[2],dhdt,basal_melting,surface_mass_balance;
     491        IssmDouble *xyz_list= NULL;
     492
     493        /*Get list of cost functions*/
     494        int *responses = NULL;
     495        int  num_responses,resp,solution;
     496        element->FindParam(&num_responses,InversionNumCostFunctionsEnum);
     497        element->FindParam(&responses,NULL,InversionCostFunctionsEnum);
     498        element->FindParam(&solution,SolutionTypeEnum);
     499        if(solution!=BalancethicknessSoftSolutionEnum) _error_("not implemented yet");
     500        if(control_type!=ThicknessEnum)                _error_("Control "<<EnumToStringx(control_type)<<" not supported");
     501
     502        /*Fetch number of vertices for this finite element*/
     503        int numvertices = element->GetNumberOfVertices();
     504
     505        /*Initialize some vectors*/
     506        IssmDouble* basis         = xNew<IssmDouble>(numvertices);
     507        IssmDouble* dbasis        = xNew<IssmDouble>(2*numvertices);
     508        IssmDouble* ge            = xNewZeroInit<IssmDouble>(numvertices);
     509        int*        vertexpidlist = xNew<int>(numvertices);
     510
     511        /*Retrieve all inputs we will be needing: */
     512        element->GetVerticesCoordinates(&xyz_list);
     513        element->GradientIndexing(&vertexpidlist[0],control_index);
     514        Input* thickness_input            = element->GetInput(ThicknessEnum);                          _assert_(thickness_input);
     515        Input* thicknessobs_input         = element->GetInput(InversionThicknessObsEnum);              _assert_(thicknessobs_input);
     516        Input* weights_input              = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     517        Input* vx_input                   = element->GetInput(VxEnum);                                 _assert_(vx_input);
     518        Input* vy_input                   = element->GetInput(VyEnum);                                 _assert_(vy_input);
     519        Input* surface_mass_balance_input = element->GetInput(SurfaceforcingsMassBalanceEnum);         _assert_(surface_mass_balance_input);
     520        Input* basal_melting_input        = element->GetInput(BasalforcingsMeltingRateEnum);           _assert_(basal_melting_input);
     521        Input* dhdt_input                 = element->GetInput(BalancethicknessThickeningRateEnum);     _assert_(dhdt_input);
     522
     523        /* Start  looping on the number of gaussian points: */
     524        Gauss* gauss=element->NewGauss(4);
     525        for(int ig=gauss->begin();ig<gauss->end();ig++){
     526                gauss->GaussPoint(ig);
     527
     528                thickness_input->GetInputValue(&thickness, gauss);
     529                thickness_input->GetInputDerivativeValue(&dH[0],xyz_list,gauss);
     530                thicknessobs_input->GetInputValue(&thicknessobs, gauss);
     531
     532                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     533                element->NodalFunctionsP1(basis,gauss);
     534                element->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss);
     535
     536                /*Deal with first part (partial derivative a J with respect to k)*/
     537                for(resp=0;resp<num_responses;resp++){
     538
     539                        weights_input->GetInputValue(&weight,gauss,responses[resp]);
     540
     541                        switch(responses[resp]){
     542                                case ThicknessAbsMisfitEnum:
     543                                        for(int i=0;i<numvertices;i++) ge[i]+= (thicknessobs-thickness)*weight*Jdet*gauss->weight*basis[i];
     544                                        break;
     545                                case ThicknessAbsGradientEnum:
     546                                        for(int i=0;i<numvertices;i++) ge[i]+= - weight*dH[0]*dbasis[0*numvertices+i]*Jdet*gauss->weight;
     547                                        for(int i=0;i<numvertices;i++) ge[i]+= - weight*dH[1]*dbasis[1*numvertices+i]*Jdet*gauss->weight;
     548                                        break;
     549                                case ThicknessAlongGradientEnum:
     550                                        vx_input->GetInputValue(&vx,gauss);
     551                                        vy_input->GetInputValue(&vy,gauss);
     552                                        vel = sqrt(vx*vx+vy*vy);
     553                                        vx  = vx/(vel+1.e-9);
     554                                        vy  = vy/(vel+1.e-9);
     555                                        for(int i=0;i<numvertices;i++) ge[i]+= - weight*(dH[0]*vx+dH[1]*vy)*(dbasis[0*numvertices+i]*vx+dbasis[1*numvertices+i]*vy)*Jdet*gauss->weight;
     556                                        break;
     557                                case ThicknessAcrossGradientEnum:
     558                                        vx_input->GetInputValue(&vx,gauss);
     559                                        vy_input->GetInputValue(&vy,gauss);
     560                                        vel = sqrt(vx*vx+vy*vy);
     561                                        vx  = vx/(vel+1.e-9);
     562                                        vy  = vy/(vel+1.e-9);
     563                                        for(int i=0;i<numvertices;i++) ge[i]+= - weight*(dH[0]*(-vy)+dH[1]*vx)*(dbasis[0*numvertices+i]*(-vy)+dbasis[1*numvertices+i]*vx)*Jdet*gauss->weight;
     564                                        break;
     565                                case BalancethicknessMisfitEnum:
     566                                        surface_mass_balance_input->GetInputValue(&surface_mass_balance,gauss);
     567                                        basal_melting_input->GetInputValue(&basal_melting,gauss);
     568                                        dhdt_input->GetInputValue(&dhdt,gauss);
     569                                        vx_input->GetInputValue(&vx,gauss);
     570                                        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
     571                                        vy_input->GetInputValue(&vy,gauss);
     572                                        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
     573                                        for(int i=0;i<numvertices;i++){
     574                                                ge[i]+= - weight*Jdet*gauss->weight*(
     575                                                        (vx*dH[0]+vy*dH[1] + thickness*(dvx[0]+dvy[1]))*(vx*dbasis[0*numvertices+i]+ vy*dbasis[1*numvertices+i] + basis[i]*(dvx[0]+dvy[1]))
     576                                                        -(surface_mass_balance-basal_melting-dhdt)*(vx*dbasis[0*numvertices+i]+ vy*dbasis[1*numvertices+i] + basis[i]*(dvx[0]+dvy[1]))
     577                                                        );
     578                                        }
     579                                        break;
     580                                default:
     581                                        _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
     582                        }
     583                }
     584        }
     585        gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
     586
     587        /*Clean up and return*/
     588        xDelete<IssmDouble>(xyz_list);
     589        xDelete<IssmDouble>(basis);
     590        xDelete<IssmDouble>(ge);
     591        xDelete<int>(vertexpidlist);
     592        xDelete<int>(responses);
     593        delete gauss;
     594
     595
    483596}/*}}}*/
    484597void BalancethicknessAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/

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

    r18060 r18062  
    277277                #endif
    278278
    279                 virtual void   Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index)=0;
    280279                virtual void   ControlInputGetGradient(Vector<IssmDouble>* gradient,int enum_type,int control_index)=0;
    281280                virtual void   ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index)=0;

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

    r18060 r18062  
    28952895
    28962896}/*}}}*/
    2897 void       Penta::Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index){/*{{{*/
    2898         /*dJ/dalpha = ∂L/∂alpha = ∂J/∂alpha + ∂/∂alpha(KU-F)*/
    2899 
    2900         int   approximation;
    2901         Tria* tria=NULL;
    2902 
    2903         /*If on water, skip grad (=0): */
    2904         if(!IsIceInElement())return;
    2905                                        
    2906         /*First deal with ∂/∂alpha(KU-F)*/
    2907         switch(control_type){
    2908 
    2909                 case FrictionCoefficientEnum:
    2910                         inputs->GetInputValue(&approximation,ApproximationEnum);
    2911                         switch(approximation){
    2912                                 case SSAApproximationEnum:
    2913                                         GradjDragSSA(gradient,control_index);
    2914                                         break;
    2915                                 case HOApproximationEnum:
    2916                                         GradjDragHO(gradient,control_index);
    2917                                         break;
    2918                                 case FSApproximationEnum:
    2919                                         GradjDragFS(gradient,control_index);
    2920                                         break;
    2921                                 case NoneApproximationEnum:
    2922                                         /*Gradient is 0*/
    2923                                         break;
    2924                                 default:
    2925                                         _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
    2926                         }
    2927                         break;
    2928 
    2929                 case MaterialsRheologyBbarEnum:
    2930                         inputs->GetInputValue(&approximation,ApproximationEnum);
    2931                         switch(approximation){
    2932                                 case SSAApproximationEnum:
    2933                                         GradjBbarSSA(gradient,control_index);
    2934                                         break;
    2935                                 case HOApproximationEnum:
    2936                                         GradjBbarHO(gradient,control_index);
    2937                                         break;
    2938                                 case FSApproximationEnum:
    2939                                         GradjBbarFS(gradient,control_index);
    2940                                         break;
    2941                                 case NoneApproximationEnum:
    2942                                         /*Gradient is 0*/
    2943                                         break;
    2944                                 default:
    2945                                         _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
    2946                         }
    2947                         break;
    2948 
    2949                 default:
    2950                         _error_("control type " << EnumToStringx(control_type) << " not supported yet: ");
    2951         }
    2952 
    2953         /*Now deal with ∂J/∂alpha*/
    2954         int *responses = NULL;
    2955         int num_responses,resp;
    2956         this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    2957         this->parameters->FindParam(&responses,NULL,InversionCostFunctionsEnum);
    2958 
    2959         for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
    2960 
    2961                 case ThicknessAbsMisfitEnum:
    2962                 case SurfaceAbsVelMisfitEnum:
    2963                 case SurfaceRelVelMisfitEnum:
    2964                 case SurfaceLogVelMisfitEnum:
    2965                 case SurfaceLogVxVyMisfitEnum:
    2966                 case SurfaceAverageVelMisfitEnum:
    2967                         /*Nothing, J does not depends on the parameter being inverted for*/
    2968                         break;
    2969                 case DragCoefficientAbsGradientEnum:
    2970                         if(IsOnBase()){
    2971                                 tria=(Tria*)SpawnTria(0,1,2);
    2972                                 tria->GradjDragGradient(gradient,control_index);
    2973                                 delete tria->material; delete tria;
    2974                         }
    2975                         break;
    2976                 case RheologyBbarAbsGradientEnum:
    2977                         if(IsOnBase()){
    2978                                 tria=(Tria*)SpawnTria(0,1,2);
    2979                                 tria->GradjBGradient(gradient,control_index);
    2980                                 delete tria->material; delete tria;
    2981                         }
    2982                         break;
    2983                 default:
    2984                         _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
    2985         }
    2986         xDelete<int>(responses);
    2987 }
    2988 /*}}}*/
    2989 void       Penta::GradjDragSSA(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    2990 
    2991         /*Gradient is 0 if on shelf or not on bed*/
    2992         if(IsFloating() || !IsOnBase()) return;
    2993 
    2994         /*Spawn tria*/
    2995         Tria* tria=(Tria*)SpawnTria(0,1,2);
    2996         tria->GradjDragSSA(gradient,control_index);
    2997         delete tria->material; delete tria;
    2998 
    2999 } /*}}}*/
    3000 void       Penta::GradjDragHO(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    3001 
    3002         printf("-------------- file: Penta.cpp line: %i\n",__LINE__);
    3003         int        i,j;
    3004         int        analysis_type;
    3005         int        vertexpidlist[NUMVERTICES];
    3006         IssmDouble vx,vy,lambda,mu,alpha_complement,Jdet;
    3007         IssmDouble bed,thickness,Neff,drag;
    3008         IssmDouble xyz_list[NUMVERTICES][3];
    3009         IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0};
    3010         IssmDouble dk[NDOF3];
    3011         IssmDouble grade_g[NUMVERTICES]={0.0};
    3012         IssmDouble grade_g_gaussian[NUMVERTICES];
    3013         IssmDouble basis[6];
    3014         Friction*  friction=NULL;
    3015         GaussPenta *gauss=NULL;
    3016 
    3017         /*Gradient is 0 if on shelf or not on bed*/
    3018         if(IsFloating() || !IsOnBase()) return;
    3019 
    3020         /*Retrieve all inputs and parameters*/
    3021         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    3022         GradientIndexing(&vertexpidlist[0],control_index);
    3023         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3024         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    3025         Input* adjointx_input=inputs->GetInput(AdjointxEnum);               _assert_(adjointx_input);
    3026         Input* adjointy_input=inputs->GetInput(AdjointyEnum);               _assert_(adjointy_input);
    3027         Input* vx_input=inputs->GetInput(VxEnum);                           _assert_(vx_input);
    3028         Input* vy_input=inputs->GetInput(VyEnum);                           _assert_(vy_input);
    3029         Input* dragcoefficient_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(dragcoefficient_input);
    3030 
    3031         /*Build frictoin element, needed later: */
    3032         friction=new Friction(this,2);
    3033 
    3034         /* Start  looping on the number of gaussian points: */
    3035         gauss=new GaussPenta(0,1,2,4);
    3036         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3037 
    3038                 gauss->GaussPoint(ig);
    3039 
    3040                 GetTriaJacobianDeterminant(&Jdet, &xyz_list_tria[0][0],gauss);
    3041                 GetNodalFunctionsP1(basis,gauss);
    3042 
    3043                 /*Build alpha_complement_list: */
    3044                 friction->GetAlphaComplement(&alpha_complement,gauss);
    3045 
    3046                 dragcoefficient_input->GetInputValue(&drag, gauss);
    3047                 adjointx_input->GetInputValue(&lambda, gauss);
    3048                 adjointy_input->GetInputValue(&mu, gauss);
    3049                 vx_input->GetInputValue(&vx,gauss);
    3050                 vy_input->GetInputValue(&vy,gauss);
    3051                 dragcoefficient_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
    3052 
    3053                 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    3054                 for (i=0;i<NUMVERTICES;i++){
    3055                         grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
    3056                 }
    3057 
    3058                 /*Add gradje_g_gaussian vector to gradje_g: */
    3059                 for(i=0;i<NUMVERTICES;i++){
    3060                         _assert_(!xIsNan<IssmDouble>(grade_g[i]));
    3061                         grade_g[i]+=grade_g_gaussian[i];
    3062                 }
    3063         }
    3064         gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL);
    3065 
    3066         /*Clean up and return*/
    3067         delete gauss;
    3068         delete friction;
    3069 }
    3070 /*}}}*/
    3071 void       Penta::GradjDragFS(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    3072 
    3073         int        i,j;
    3074         int        analysis_type;
    3075         int        vertexpidlist[NUMVERTICES];
    3076         IssmDouble bed,thickness,Neff;
    3077         IssmDouble lambda,mu,xi,Jdet,vx,vy,vz;
    3078         IssmDouble alpha_complement,drag;
    3079         IssmDouble surface_normal[3],bed_normal[3];
    3080         IssmDouble xyz_list[NUMVERTICES][3];
    3081         IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0};
    3082         IssmDouble dk[NDOF3];
    3083         IssmDouble basis[6];
    3084         IssmDouble grade_g[NUMVERTICES]={0.0};
    3085         IssmDouble grade_g_gaussian[NUMVERTICES];
    3086         Friction*  friction=NULL;
    3087         GaussPenta* gauss=NULL;
    3088 
    3089         /*Gradient is 0 if on shelf or not on bed*/
    3090         if(IsFloating() || !IsOnBase()) return;
    3091 
    3092         /*Retrieve all inputs and parameters*/
    3093         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    3094         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3095         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    3096         GradientIndexing(&vertexpidlist[0],control_index);
    3097         Input* drag_input    =inputs->GetInput(FrictionCoefficientEnum); _assert_(drag_input);
    3098         Input* vx_input      =inputs->GetInput(VxEnum);                  _assert_(vx_input);
    3099         Input* vy_input      =inputs->GetInput(VyEnum);                  _assert_(vy_input);
    3100         Input* vz_input      =inputs->GetInput(VzEnum);                  _assert_(vz_input);
    3101         Input* adjointx_input=inputs->GetInput(AdjointxEnum);            _assert_(adjointx_input);
    3102         Input* adjointy_input=inputs->GetInput(AdjointyEnum);            _assert_(adjointy_input);
    3103         Input* adjointz_input=inputs->GetInput(AdjointzEnum);            _assert_(adjointz_input);
    3104 
    3105         /*Build frictoin element, needed later: */
    3106         friction=new Friction(this,3);
    3107 
    3108         /* Start  looping on the number of gaussian points: */
    3109         gauss=new GaussPenta(0,1,2,4);
    3110         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3111 
    3112                 gauss->GaussPoint(ig);
    3113 
    3114                 /*Recover alpha_complement and drag: */
    3115                 friction->GetAlphaComplement(&alpha_complement,gauss);
    3116                 drag_input->GetInputValue(&drag,gauss);
    3117 
    3118                 /*recover lambda mu and xi: */
    3119                 adjointx_input->GetInputValue(&lambda,gauss);
    3120                 adjointy_input->GetInputValue(&mu    ,gauss);
    3121                 adjointz_input->GetInputValue(&xi    ,gauss);
    3122 
    3123                 /*recover vx vy and vz: */
    3124                 vx_input->GetInputValue(&vx, gauss);
    3125                 vy_input->GetInputValue(&vy, gauss);
    3126                 vz_input->GetInputValue(&vz, gauss);
    3127 
    3128                 /*Get normal vector to the bed */
    3129                 NormalTop(&surface_normal[0],&xyz_list_tria[0][0]);
    3130 
    3131                 bed_normal[0]=-surface_normal[0]; //Function is for upper surface, so the normal to the bed is the opposite of the result
    3132                 bed_normal[1]=-surface_normal[1];
    3133                 bed_normal[2]=-surface_normal[2];
    3134 
    3135                 /* Get Jacobian determinant: */
    3136                 GetTriaJacobianDeterminant(&Jdet,&xyz_list_tria[0][0],gauss);
    3137                 GetNodalFunctionsP1(basis, gauss);
    3138 
    3139                 /*Get k derivative: dk/dx */
    3140                 drag_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
    3141 
    3142                 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    3143                 for (i=0;i<NUMVERTICES;i++){
    3144                         //standard gradient dJ/dki
    3145                         grade_g_gaussian[i]=(
    3146                                                 -lambda*(2*drag*alpha_complement*(vx - vz*bed_normal[0]*bed_normal[2]))
    3147                                                 -mu    *(2*drag*alpha_complement*(vy - vz*bed_normal[1]*bed_normal[2]))
    3148                                                 -xi    *(2*drag*alpha_complement*(-vx*bed_normal[0]*bed_normal[2]-vy*bed_normal[1]*bed_normal[2]))
    3149                                                 )*Jdet*gauss->weight*basis[i];
    3150                 }
    3151 
    3152                 /*Add gradje_g_gaussian vector to gradje_g: */
    3153                 for( i=0; i<NUMVERTICES; i++)grade_g[i]+=grade_g_gaussian[i];
    3154         }
    3155 
    3156         gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL);
    3157 
    3158         delete friction;
    3159         delete gauss;
    3160 }
    3161 /*}}}*/
    3162 void       Penta::GradjBbarSSA(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    3163 
    3164         /*This element should be collapsed into a tria element at its base*/
    3165         if (!IsOnBase()) return;
    3166 
    3167         /*Depth Average B*/
    3168         this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
    3169         if(this->material->IsDamage())this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum);
    3170 
    3171         /*Collapse element to the base*/
    3172         Tria* tria=(Tria*)SpawnTria(0,1,2);
    3173         tria->GradjBSSA(gradient,control_index);
    3174         delete tria->material; delete tria;
    3175 
    3176         /*delete Average B*/
    3177         this->inputs->DeleteInput(MaterialsRheologyBbarEnum);
    3178         this->inputs->DeleteInput(DamageDbarEnum);
    3179 
    3180 } /*}}}*/
    3181 void       Penta::GradjBbarHO(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    3182 
    3183         /*Gradient is computed on bed only (Bbar)*/
    3184         if (!IsOnBase()) return;
    3185 
    3186         /*Depth Average B and D*/
    3187         this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
    3188         if(this->material->IsDamage())this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum);
    3189 
    3190         /*Collapse element to the base*/
    3191         Tria* tria=(Tria*)SpawnTria(0,1,2);
    3192         tria->GradjBSSA(gradient,control_index);    //We use SSA as an estimate for now
    3193         delete tria->material; delete tria;
    3194 
    3195         /*delete Average B*/
    3196         this->inputs->DeleteInput(MaterialsRheologyBbarEnum);
    3197         this->inputs->DeleteInput(DamageDbarEnum);
    3198 } /*}}}*/
    3199 void       Penta::GradjBbarFS(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    3200 
    3201         /*Gradient is computed on bed only (Bbar)*/
    3202         if (!IsOnBase()) return;
    3203 
    3204         /*Depth Average B and D*/
    3205         this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
    3206         if(this->material->IsDamage())this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum);
    3207 
    3208         /*Collapse element to the base*/
    3209         Tria* tria=(Tria*)SpawnTria(0,1,2);
    3210         tria->GradjBSSA(gradient,control_index);    //We use SSA as an estimate for now
    3211         delete tria->material; delete tria;
    3212 
    3213         /*delete Average B*/
    3214         this->inputs->DeleteInput(MaterialsRheologyBbarEnum);
    3215         this->inputs->DeleteInput(DamageDbarEnum);
    3216 } /*}}}*/
    32172897void       Penta::InputControlUpdate(IssmDouble scalar,bool save_parameter){/*{{{*/
    32182898

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

    r18058 r18062  
    131131                #endif
    132132
    133                 void   Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index);
    134                 void   GradjDragSSA(Vector<IssmDouble>* gradient,int control_index);
    135                 void   GradjDragHO(Vector<IssmDouble>* gradient,int control_index);
    136                 void   GradjDragFS(Vector<IssmDouble>* gradient,int control_index);
    137                 void   GradjBbarSSA(Vector<IssmDouble>* gradient,int control_index);
    138                 void   GradjBbarHO(Vector<IssmDouble>* gradient,int control_index);
    139                 void   GradjBbarFS(Vector<IssmDouble>* gradient,int control_index);
    140133                void   GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data);
    141134                void   SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index);

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

    r18058 r18062  
    116116
    117117}/*}}}*/
    118 void       Seg::GradjDragFS(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    119 
    120         int        i;
    121         int        analysis_type;
    122         int        vertexpidlist[NUMVERTICES];
    123         int        connectivity[NUMVERTICES];
    124         IssmDouble vx,lambda,alpha_complement,drag,Jdet;
    125         IssmDouble xyz_list[NUMVERTICES][3];
    126         IssmDouble dk[NDOF2];
    127         IssmDouble grade_g[NUMVERTICES]={0.0};
    128         IssmDouble grade_g_gaussian[NUMVERTICES];
    129         IssmDouble basis[3];
    130         Friction*  friction=NULL;
    131         GaussSeg  *gauss=NULL;
    132 
    133         if(IsFloating())return;
    134 
    135         /*retrive parameters: */
    136         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    137         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    138         GradientIndexing(&vertexpidlist[0],control_index);
    139         this->GetVerticesConnectivityList(&connectivity[0]);
    140 
    141         /*Build frictoin element, needed later: */
    142         friction=new Friction(this,1);
    143 
    144         /*Retrieve all inputs we will be needing: */
    145         Input* adjointx_input=inputs->GetInput(AdjointxEnum);                   _assert_(adjointx_input);
    146         Input* vx_input=inputs->GetInput(VxEnum);                               _assert_(vx_input);
    147         Input* dragcoefficient_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(dragcoefficient_input);
    148 
    149         /* Start  looping on the number of gaussian points: */
    150         gauss=new GaussSeg(4);
    151         for(int ig=gauss->begin();ig<gauss->end();ig++){
    152 
    153                 gauss->GaussPoint(ig);
    154 
    155                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    156                 GetNodalFunctions(basis, gauss);
    157 
    158                 /*Build alpha_complement_list: */
    159                 friction->GetAlphaComplement(&alpha_complement,gauss);
    160 
    161                 dragcoefficient_input->GetInputValue(&drag, gauss);
    162                 adjointx_input->GetInputValue(&lambda, gauss);
    163                 vx_input->GetInputValue(&vx,gauss);
    164                 dragcoefficient_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
    165 
    166                 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    167                 for (i=0;i<NUMVERTICES;i++){
    168                         grade_g_gaussian[i]=-2*drag*alpha_complement*(lambda*vx)*Jdet*gauss->weight*basis[i];
    169                 }
    170 
    171                 /*Add gradje_g_gaussian vector to gradje_g: */
    172                 for(i=0;i<NUMVERTICES;i++){
    173                         _assert_(!xIsNan<IssmDouble>(grade_g[i]));
    174                         grade_g[i]+=grade_g_gaussian[i];
    175                 }
    176         }
    177         gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL);
    178 
    179         /*Clean up and return*/
    180         delete gauss;
    181         delete friction;
    182 }
    183 /*}}}*/
    184 void       Seg::GradjDragGradient(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    185 
    186         int        i;
    187         int        vertexpidlist[NUMVERTICES];
    188         IssmDouble Jdet,weight;
    189         IssmDouble xyz_list[NUMVERTICES][3];
    190         IssmDouble dbasis[NDOF2][NUMVERTICES];
    191         IssmDouble dk[NDOF2];
    192         IssmDouble grade_g[NUMVERTICES]={0.0};
    193         GaussSeg  *gauss=NULL;
    194 
    195         /*Retrieve all inputs we will be needing: */
    196         if(IsFloating())return;
    197         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    198         GradientIndexing(&vertexpidlist[0],control_index);
    199         Input* dragcoefficient_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(dragcoefficient_input);
    200         Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);                 _assert_(weights_input);
    201 
    202         /* Start  looping on the number of gaussian points: */
    203         gauss=new GaussSeg(2);
    204         for(int ig=gauss->begin();ig<gauss->end();ig++){
    205 
    206                 gauss->GaussPoint(ig);
    207 
    208                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    209                 GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
    210                 weights_input->GetInputValue(&weight,gauss,DragCoefficientAbsGradientEnum);
    211 
    212                 /*Build alpha_complement_list: */
    213                 dragcoefficient_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
    214 
    215                 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    216                 for (i=0;i<NUMVERTICES;i++){
    217                         grade_g[i]+=-weight*Jdet*gauss->weight*dbasis[0][i]*dk[0];
    218                         _assert_(!xIsNan<IssmDouble>(grade_g[i]));
    219                 }
    220         }
    221         gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL);
    222 
    223         /*Clean up and return*/
    224         delete gauss;
    225 }
    226 /*}}}*/
    227118bool       Seg::IsIcefront(void){/*{{{*/
    228119

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

    r18058 r18062  
    167167#endif
    168168
    169                 void       Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index){_error_("not implemented yet");};
    170                 void       GradjBGradient(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
    171                 void       GradjDGradient(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
    172                 void       GradjBSSA(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
    173                 void       GradjDSSA(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
    174                 void       GradjDragSSA(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
    175                 void       GradjDragFS(Vector<IssmDouble>* gradient,int control_index);
    176                 void       GradjDragGradient(Vector<IssmDouble>* gradient,int control_index);
    177                 void       GradjDhDtBalancedthickness(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
    178                 void       GradjVxBalancedthickness(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
    179                 void       GradjVyBalancedthickness(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
    180                 void       GradjThicknessBalancethicknessSoft(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
    181169                void       GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data){_error_("not implemented yet");};
    182170                void       SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){_error_("not implemented yet");};

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

    r18001 r18062  
    173173                IssmDouble DragCoefficientAbsGradient(void){_error_("not implemented yet");};
    174174                void       GradientIndexing(int* indexing,int control_index){_error_("not implemented yet");};
    175                 void       Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index){_error_("not implemented yet");};
    176                 void       GradjBGradient(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
    177                 void       GradjDGradient(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
    178                 void       GradjBSSA(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
    179                 void       GradjDSSA(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
    180                 void       GradjDragSSA(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
    181                 void       GradjDragFS(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
    182                 void       GradjDragGradient(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
    183                 void       GradjDhDtBalancedthickness(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
    184                 void       GradjVxBalancedthickness(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
    185                 void       GradjVyBalancedthickness(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
    186                 void       GradjThicknessBalancethicknessSoft(Vector<IssmDouble>* gradient,int control_index){_error_("not implemented yet");};
    187175                void       GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data){_error_("not implemented yet");};
    188176                void       SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){_error_("not implemented yet");};

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

    r18060 r18062  
    27732773
    27742774}/*}}}*/
    2775 void       Tria::Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index){/*{{{*/
    2776         /*dJ/dalpha = ∂L/∂alpha = ∂J/∂alpha + ∂/∂alpha(KU-F)*/
    2777 
    2778         int   approximation;
    2779         Seg*  seg=NULL;
    2780 
    2781         /*If on water, grad = 0: */
    2782         if(!IsIceInElement()) return;
    2783 
    2784         /*First deal with ∂/∂alpha(KU-F)*/
    2785         switch(control_type){
    2786                 case FrictionCoefficientEnum:
    2787                         inputs->GetInputValue(&approximation,ApproximationEnum);
    2788                         switch(approximation){
    2789                                 case SSAApproximationEnum:
    2790                                         GradjDragSSA(gradient,control_index);
    2791                                         break;
    2792                                 case FSApproximationEnum:
    2793                                         GradjDragFS(gradient,control_index);
    2794                                         break;
    2795                                 case NoneApproximationEnum:
    2796                                         /*Gradient is 0*/
    2797                                         break;
    2798                                 default:
    2799                                         _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
    2800                         }
    2801                         break;
    2802                 case MaterialsRheologyBbarEnum:
    2803                         GradjBSSA(gradient,control_index);
    2804                         break;
    2805                 case DamageDbarEnum:
    2806                         GradjDSSA(gradient,control_index);
    2807                         break;
    2808                 case BalancethicknessThickeningRateEnum:
    2809                         GradjDhDtBalancedthickness(gradient,control_index);
    2810                         break;
    2811                 case VxEnum:
    2812                         GradjVxBalancedthickness(gradient,control_index);
    2813                         break;
    2814                 case VyEnum:
    2815                         GradjVyBalancedthickness(gradient,control_index);
    2816                         break;
    2817                 case ThicknessEnum:
    2818                         GradjThicknessBalancethicknessSoft(gradient,control_index);
    2819                         break;
    2820                 case BalancethicknessApparentMassbalanceEnum:
    2821                         GradjAdotBulancethickness2(gradient,control_index);
    2822                         break;
    2823                 default:
    2824                         _error_("control type not supported yet: " << EnumToStringx(control_type));
    2825         }
    2826 
    2827         /*Now deal with ∂J/∂alpha*/
    2828         int        *responses = NULL;
    2829         int         num_responses,resp;
    2830         this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    2831         this->parameters->FindParam(&responses,NULL,InversionCostFunctionsEnum);
    2832 
    2833         for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
    2834                 //FIXME: the control type should be checked somewhere (with respect to what variable are we taking the gradient!)
    2835 
    2836                 case ThicknessAbsMisfitEnum:
    2837                 case ThicknessAbsGradientEnum:
    2838                 case ThicknessAlongGradientEnum:
    2839                 case ThicknessAcrossGradientEnum:
    2840                 case BalancethicknessMisfitEnum:
    2841                 case Balancethickness2MisfitEnum:
    2842                 case SurfaceAbsVelMisfitEnum:
    2843                 case SurfaceRelVelMisfitEnum:
    2844                 case SurfaceLogVelMisfitEnum:
    2845                 case SurfaceLogVxVyMisfitEnum:
    2846                 case SurfaceAverageVelMisfitEnum:
    2847                         /*Nothing, J does not depends on the parameter being inverted for*/
    2848                         break;
    2849                 case DragCoefficientAbsGradientEnum:
    2850                         inputs->GetInputValue(&approximation,ApproximationEnum);
    2851                         switch(approximation){
    2852                                 case SSAApproximationEnum:
    2853                                         GradjDragGradient(gradient,control_index);
    2854                                         break;
    2855                                 case FSApproximationEnum:{
    2856                                         if(IsFloating() || !IsOnBase()) return;
    2857                                         int index1,index2;
    2858                                         this->EdgeOnBaseIndices(&index1,&index2);
    2859                                         Seg* seg = SpawnSeg(index1,index2);
    2860                                         seg->GradjDragGradient(gradient,control_index);
    2861                                         seg->DeleteMaterials(); delete seg;
    2862                                         break;
    2863                                                                                                  }
    2864                                 case NoneApproximationEnum:
    2865                                         /*Gradient is 0*/
    2866                                         break;
    2867                                 default:
    2868                                         _error_("approximation " << EnumToStringx(approximation) << " not supported yet");
    2869                         }
    2870                         break;
    2871                 case RheologyBbarAbsGradientEnum:
    2872                         GradjBGradient(gradient,control_index);
    2873                         break;
    2874                 default:
    2875                         _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
    2876         }
    2877 
    2878         xDelete<int>(responses);
    2879 }
    2880 /*}}}*/
    2881 void       Tria::GradjBGradient(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    2882 
    2883         int        i;
    2884         int        vertexpidlist[NUMVERTICES];
    2885         IssmDouble Jdet,weight;
    2886         IssmDouble xyz_list[NUMVERTICES][3];
    2887         IssmDouble dbasis[NDOF2][NUMVERTICES];
    2888         IssmDouble dk[NDOF2];
    2889         IssmDouble grade_g[NUMVERTICES]={0.0};
    2890         GaussTria  *gauss=NULL;
    2891 
    2892         /*Retrieve all inputs we will be needing: */
    2893         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    2894         GradientIndexing(&vertexpidlist[0],control_index);
    2895         Input* rheologyb_input=inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
    2896         Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);                _assert_(weights_input);
    2897 
    2898         /* Start  looping on the number of gaussian points: */
    2899         gauss=new GaussTria(2);
    2900         for(int ig=gauss->begin();ig<gauss->end();ig++){
    2901 
    2902                 gauss->GaussPoint(ig);
    2903 
    2904                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    2905                 GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
    2906                 weights_input->GetInputValue(&weight,gauss,RheologyBbarAbsGradientEnum);
    2907 
    2908                 /*Build alpha_complement_list: */
    2909                 rheologyb_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
    2910 
    2911                 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    2912                 for (i=0;i<NUMVERTICES;i++) grade_g[i]+=-weight*Jdet*gauss->weight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]);
    2913         }
    2914         gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL);
    2915 
    2916         /*Clean up and return*/
    2917         delete gauss;
    2918 }
    2919 /*}}}*/
    2920 void       Tria::GradjBSSA(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    2921 
    2922         /*Intermediaries*/
    2923         int        i;
    2924         int        doflist[NUMVERTICES];
    2925         IssmDouble vx,vy,lambda,mu,thickness,Jdet;
    2926         IssmDouble viscosity_complement;
    2927         IssmDouble dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dB[NDOF2];
    2928         IssmDouble xyz_list[NUMVERTICES][3];
    2929         IssmDouble basis[3],epsilon[3];
    2930         IssmDouble grad[NUMVERTICES]={0.0};
    2931         GaussTria *gauss = NULL;
    2932 
    2933         /* Get node coordinates and dof list: */
    2934         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    2935         GradientIndexing(&doflist[0],control_index);
    2936 
    2937         /*Retrieve all inputs*/
    2938         Input* thickness_input=inputs->GetInput(ThicknessEnum);                     _assert_(thickness_input);
    2939         Input* vx_input=inputs->GetInput(VxEnum);                                   _assert_(vx_input);
    2940         Input* vy_input=inputs->GetInput(VyEnum);                                   _assert_(vy_input);
    2941         Input* adjointx_input=inputs->GetInput(AdjointxEnum);                       _assert_(adjointx_input);
    2942         Input* adjointy_input=inputs->GetInput(AdjointyEnum);                       _assert_(adjointy_input);
    2943         Input* rheologyb_input=inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
    2944 
    2945         /* Start  looping on the number of gaussian points: */
    2946         gauss=new GaussTria(4);
    2947         for(int ig=gauss->begin();ig<gauss->end();ig++){
    2948 
    2949                 gauss->GaussPoint(ig);
    2950 
    2951                 thickness_input->GetInputValue(&thickness,gauss);
    2952                 rheologyb_input->GetInputDerivativeValue(&dB[0],&xyz_list[0][0],gauss);
    2953                 vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
    2954                 vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
    2955                 adjointx_input->GetInputDerivativeValue(&dadjx[0],&xyz_list[0][0],gauss);
    2956                 adjointy_input->GetInputDerivativeValue(&dadjy[0],&xyz_list[0][0],gauss);
    2957 
    2958                 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    2959                 material->GetViscosityComplement(&viscosity_complement,&epsilon[0]);
    2960 
    2961                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    2962                 GetNodalFunctions(basis,gauss);
    2963 
    2964                 /*standard gradient dJ/dki*/
    2965                 for (i=0;i<NUMVERTICES;i++) grad[i]+=-viscosity_complement*thickness*(
    2966                                         (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
    2967                                         )*Jdet*gauss->weight*basis[i];
    2968         }
    2969 
    2970         gradient->SetValues(NUMVERTICES,doflist,grad,ADD_VAL);
    2971 
    2972         /*clean-up*/
    2973         delete gauss;
    2974 }
    2975 /*}}}*/
    2976 void       Tria::GradjDSSA(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    2977 
    2978         /*Intermediaries*/
    2979         int        i;
    2980         int        doflist[NUMVERTICES];
    2981         IssmDouble vx,vy,lambda,mu,thickness,Jdet;
    2982         IssmDouble viscosity_complement;
    2983         IssmDouble dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2];
    2984         IssmDouble xyz_list[NUMVERTICES][3];
    2985         IssmDouble basis[3],epsilon[3];
    2986         IssmDouble grad[NUMVERTICES]={0.0};
    2987         GaussTria *gauss = NULL;
    2988 
    2989         /* Get node coordinates and dof list: */
    2990         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    2991         GradientIndexing(&doflist[0],control_index);
    2992 
    2993         /*Retrieve all inputs*/
    2994         Input* thickness_input=inputs->GetInput(ThicknessEnum);                     _assert_(thickness_input);
    2995         Input* vx_input=inputs->GetInput(VxEnum);                                   _assert_(vx_input);
    2996         Input* vy_input=inputs->GetInput(VyEnum);                                   _assert_(vy_input);
    2997         Input* adjointx_input=inputs->GetInput(AdjointxEnum);                       _assert_(adjointx_input);
    2998         Input* adjointy_input=inputs->GetInput(AdjointyEnum);                       _assert_(adjointy_input);
    2999         if(this->material->IsDamage()){
    3000                 Input* rheologyd_input=inputs->GetInput(DamageDbarEnum); _assert_(rheologyd_input);
    3001         }
    3002 
    3003         /* Start  looping on the number of gaussian points: */
    3004         gauss=new GaussTria(4);
    3005         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3006 
    3007                 gauss->GaussPoint(ig);
    3008 
    3009                 thickness_input->GetInputValue(&thickness,gauss);
    3010                 vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
    3011                 vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
    3012                 adjointx_input->GetInputDerivativeValue(&dadjx[0],&xyz_list[0][0],gauss);
    3013                 adjointy_input->GetInputDerivativeValue(&dadjy[0],&xyz_list[0][0],gauss);
    3014 
    3015                 this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    3016                 material->GetViscosityDComplement(&viscosity_complement,&epsilon[0]);
    3017 
    3018                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3019                 GetNodalFunctions(basis,gauss);
    3020 
    3021                 /*standard gradient dJ/dki*/
    3022                 for (i=0;i<NUMVERTICES;i++) grad[i]+=-viscosity_complement*thickness*(
    3023                                         (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
    3024                                         )*Jdet*gauss->weight*basis[i];
    3025         }
    3026 
    3027         gradient->SetValues(NUMVERTICES,doflist,grad,ADD_VAL);
    3028 
    3029         /*clean-up*/
    3030         delete gauss;
    3031 }
    3032 /*}}}*/
    3033 void       Tria::GradjDragSSA(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    3034 
    3035         int        i;
    3036         int        analysis_type;
    3037         int        vertexpidlist[NUMVERTICES];
    3038         int        connectivity[NUMVERTICES];
    3039         IssmDouble vx,vy,lambda,mu,alpha_complement,Jdet;
    3040         IssmDouble bed,thickness,Neff,drag;
    3041         IssmDouble xyz_list[NUMVERTICES][3];
    3042         IssmDouble dk[NDOF2];
    3043         IssmDouble grade_g[NUMVERTICES]={0.0};
    3044         IssmDouble grade_g_gaussian[NUMVERTICES];
    3045         IssmDouble basis[3];
    3046         Friction*  friction=NULL;
    3047         GaussTria  *gauss=NULL;
    3048 
    3049         if(IsFloating())return;
    3050 
    3051         /*retrive parameters: */
    3052         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    3053         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3054         GradientIndexing(&vertexpidlist[0],control_index);
    3055         this->GetVerticesConnectivityList(&connectivity[0]);
    3056 
    3057         /*Build frictoin element, needed later: */
    3058         friction=new Friction(this,2);
    3059 
    3060         /*Retrieve all inputs we will be needing: */
    3061         Input* adjointx_input=inputs->GetInput(AdjointxEnum);                   _assert_(adjointx_input);
    3062         Input* adjointy_input=inputs->GetInput(AdjointyEnum);                   _assert_(adjointy_input);
    3063         Input* vx_input=inputs->GetInput(VxEnum);                               _assert_(vx_input);
    3064         Input* vy_input=inputs->GetInput(VyEnum);                               _assert_(vy_input);
    3065         Input* dragcoefficient_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(dragcoefficient_input);
    3066 
    3067         /* Start  looping on the number of gaussian points: */
    3068         gauss=new GaussTria(4);
    3069         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3070 
    3071                 gauss->GaussPoint(ig);
    3072 
    3073                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3074                 GetNodalFunctions(basis, gauss);
    3075 
    3076                 /*Build alpha_complement_list: */
    3077                 friction->GetAlphaComplement(&alpha_complement,gauss);
    3078 
    3079                 dragcoefficient_input->GetInputValue(&drag, gauss);
    3080                 adjointx_input->GetInputValue(&lambda, gauss);
    3081                 adjointy_input->GetInputValue(&mu, gauss);
    3082                 vx_input->GetInputValue(&vx,gauss);
    3083                 vy_input->GetInputValue(&vy,gauss);
    3084                 dragcoefficient_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
    3085 
    3086                 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    3087                 for (i=0;i<NUMVERTICES;i++){
    3088                         grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
    3089                 }
    3090 
    3091                 /*Add gradje_g_gaussian vector to gradje_g: */
    3092                 for(i=0;i<NUMVERTICES;i++){
    3093                         _assert_(!xIsNan<IssmDouble>(grade_g[i]));
    3094                         grade_g[i]+=grade_g_gaussian[i];
    3095                 }
    3096         }
    3097         /*Analytical gradient*/
    3098         //delete gauss;
    3099         //gauss=new GaussTria();
    3100         //for (int iv=0;iv<NUMVERTICES;iv++){
    3101         //      gauss->GaussVertex(iv);
    3102         //      friction->GetAlphaComplement(&alpha_complement,gauss);
    3103         //      dragcoefficient_input->GetInputValue(&drag, gauss);
    3104         //      adjointx_input->GetInputValue(&lambda, gauss);
    3105         //      adjointy_input->GetInputValue(&mu, gauss);
    3106         //      vx_input->GetInputValue(&vx,gauss);
    3107         //      vy_input->GetInputValue(&vy,gauss);
    3108         //      grade_g[iv] = -2*1.e+7*drag*alpha_complement*(lambda*vx+mu*vy)/((IssmDouble)connectivity[iv]);
    3109         //}
    3110         /*End Analytical gradient*/
    3111 
    3112         gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL);
    3113 
    3114         /*Clean up and return*/
    3115         delete gauss;
    3116         delete friction;
    3117 }
    3118 /*}}}*/
    3119 void       Tria::GradjDragFS(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    3120 
    3121         /*Gradient is 0 if on shelf or not on bed*/
    3122         if(IsFloating() || !IsOnBase()) return;
    3123 
    3124         int index1,index2;
    3125         this->EdgeOnBaseIndices(&index1,&index2);
    3126         Seg* seg = SpawnSeg(index1,index2);
    3127         seg->GradjDragFS(gradient,control_index);
    3128         seg->DeleteMaterials(); delete seg;
    3129 }
    3130 /*}}}*/
    3131 void       Tria::GradjDragGradient(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    3132 
    3133         int        i;
    3134         int        vertexpidlist[NUMVERTICES];
    3135         IssmDouble Jdet,weight;
    3136         IssmDouble xyz_list[NUMVERTICES][3];
    3137         IssmDouble dbasis[NDOF2][NUMVERTICES];
    3138         IssmDouble dk[NDOF2];
    3139         IssmDouble grade_g[NUMVERTICES]={0.0};
    3140         GaussTria  *gauss=NULL;
    3141 
    3142         /*Retrieve all inputs we will be needing: */
    3143         if(IsFloating())return;
    3144         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3145         GradientIndexing(&vertexpidlist[0],control_index);
    3146         Input* dragcoefficient_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(dragcoefficient_input);
    3147         Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);                 _assert_(weights_input);
    3148 
    3149         /* Start  looping on the number of gaussian points: */
    3150         gauss=new GaussTria(2);
    3151         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3152 
    3153                 gauss->GaussPoint(ig);
    3154 
    3155                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3156                 GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
    3157                 weights_input->GetInputValue(&weight,gauss,DragCoefficientAbsGradientEnum);
    3158 
    3159                 /*Build alpha_complement_list: */
    3160                 dragcoefficient_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
    3161 
    3162                 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    3163                 for (i=0;i<NUMVERTICES;i++){
    3164                         grade_g[i]+=-weight*Jdet*gauss->weight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]);
    3165                         _assert_(!xIsNan<IssmDouble>(grade_g[i]));
    3166                 }
    3167         }
    3168         gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL);
    3169 
    3170         /*Clean up and return*/
    3171         delete gauss;
    3172 }
    3173 /*}}}*/
    3174 void       Tria::GradjDhDtBalancedthickness(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    3175 
    3176         /*Intermediaries*/
    3177         int    vertexpidlist[NUMVERTICES];
    3178         IssmDouble lambda[NUMVERTICES];
    3179         IssmDouble gradient_g[NUMVERTICES];
    3180 
    3181         /*Compute Gradient*/
    3182         GradientIndexing(&vertexpidlist[0],control_index);
    3183         GetInputListOnVertices(&lambda[0],AdjointEnum);
    3184         for(int i=0;i<NUMVERTICES;i++) gradient_g[i]=-lambda[i];
    3185 
    3186         gradient->SetValues(NUMVERTICES,vertexpidlist,gradient_g,INS_VAL);
    3187 }
    3188 /*}}}*/
    3189 void       Tria::GradjAdotBulancethickness2(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    3190 
    3191         /*Intermediaries*/
    3192         IssmDouble lambda,potential,weight,Jdet;
    3193         IssmDouble grade_g[NUMVERTICES];
    3194         IssmDouble basis[NUMVERTICES];
    3195         IssmDouble xyz_list[NUMVERTICES][3];
    3196         int        vertexpidlist[NUMVERTICES];
    3197 
    3198         /*Recover inputs*/
    3199         Input* lambda_input  = inputs->GetInput(AdjointEnum);                            _assert_(lambda_input);
    3200         Input* weights_input = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    3201 
    3202         /* Get node coordinates and dof list: */
    3203         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3204         GradientIndexing(&vertexpidlist[0],control_index);
    3205 
    3206         /* Start  looping on the number of gaussian points: */
    3207         Gauss* gauss=new GaussTria(2);
    3208         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3209 
    3210                 gauss->GaussPoint(ig);
    3211 
    3212                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3213                 GetNodalFunctions(&basis[0],gauss);
    3214                 weights_input->GetInputValue(&weight,gauss,Balancethickness2MisfitEnum);
    3215                 lambda_input->GetInputValue(&lambda,gauss);
    3216 
    3217                 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    3218                 for(int i=0;i<NUMVERTICES;i++){
    3219                         grade_g[i]+= - weight*Jdet*gauss->weight*basis[i]*lambda;
    3220                 }
    3221         }
    3222         gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL);
    3223 
    3224         /*Clean up and return*/
    3225         delete gauss;
    3226 }
    3227 /*}}}*/
    3228 void       Tria::GradjVxBalancedthickness(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    3229 
    3230         /*Intermediaries*/
    3231         int        i;
    3232         int        vertexpidlist[NUMVERTICES];
    3233         IssmDouble thickness,Jdet;
    3234         IssmDouble basis[3];
    3235         IssmDouble Dlambda[2],dp[2];
    3236         IssmDouble xyz_list[NUMVERTICES][3];
    3237         IssmDouble grade_g[NUMVERTICES] = {0.0};
    3238         GaussTria *gauss                = NULL;
    3239 
    3240         /* Get node coordinates and dof list: */
    3241         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3242         GradientIndexing(&vertexpidlist[0],control_index);
    3243 
    3244         /*Retrieve all inputs we will be needing: */
    3245         Input* adjoint_input=inputs->GetInput(AdjointEnum);     _assert_(adjoint_input);
    3246         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    3247 
    3248         /* Start  looping on the number of gaussian points: */
    3249         gauss=new GaussTria(2);
    3250         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3251 
    3252                 gauss->GaussPoint(ig);
    3253 
    3254                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3255                 GetNodalFunctions(basis, gauss);
    3256 
    3257                 adjoint_input->GetInputDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss);
    3258                 thickness_input->GetInputValue(&thickness, gauss);
    3259                 thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
    3260 
    3261                 for(i=0;i<NUMVERTICES;i++) grade_g[i]+=thickness*Dlambda[0]*Jdet*gauss->weight*basis[i];
    3262         }
    3263 
    3264         gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL);
    3265 
    3266         /*Clean up and return*/
    3267         delete gauss;
    3268 }
    3269 /*}}}*/
    3270 void       Tria::GradjVyBalancedthickness(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    3271 
    3272         /*Intermediaries*/
    3273         int        i;
    3274         int        vertexpidlist[NUMVERTICES];
    3275         IssmDouble thickness,Jdet;
    3276         IssmDouble basis[3];
    3277         IssmDouble Dlambda[2],dp[2];
    3278         IssmDouble xyz_list[NUMVERTICES][3];
    3279         IssmDouble grade_g[NUMVERTICES] = {0.0};
    3280         GaussTria *gauss                = NULL;
    3281 
    3282         /* Get node coordinates and dof list: */
    3283         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3284         GradientIndexing(&vertexpidlist[0],control_index);
    3285 
    3286         /*Retrieve all inputs we will be needing: */
    3287         Input* adjoint_input=inputs->GetInput(AdjointEnum);     _assert_(adjoint_input);
    3288         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    3289 
    3290         /* Start  looping on the number of gaussian points: */
    3291         gauss=new GaussTria(2);
    3292         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3293 
    3294                 gauss->GaussPoint(ig);
    3295 
    3296                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3297                 GetNodalFunctions(basis, gauss);
    3298 
    3299                 adjoint_input->GetInputDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss);
    3300                 thickness_input->GetInputValue(&thickness, gauss);
    3301                 thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
    3302 
    3303                 for(i=0;i<NUMVERTICES;i++) grade_g[i]+=thickness*Dlambda[1]*Jdet*gauss->weight*basis[i];
    3304         }
    3305         gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL);
    3306 
    3307         /*Clean up and return*/
    3308         delete gauss;
    3309 }
    3310 /*}}}*/
    3311 void       Tria::GradjThicknessBalancethicknessSoft(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    3312 
    3313         /*Intermediaries */
    3314         int         i,resp;
    3315         int         vertexpidlist[NUMVERTICES];
    3316         IssmDouble  Jdet;
    3317         IssmDouble  thickness,thicknessobs,weight;
    3318         int         num_responses;
    3319         IssmDouble  xyz_list[NUMVERTICES][3];
    3320         IssmDouble  basis[3];
    3321         IssmDouble  dbasis[NDOF2][NUMVERTICES];
    3322         IssmDouble  dH[2];
    3323         IssmDouble  vx,vy,vel;
    3324         IssmDouble  dvx[2],dvy[2];
    3325         IssmDouble dhdt,basal_melting,surface_mass_balance;
    3326         GaussTria *gauss     = NULL;
    3327         int       *responses = NULL;
    3328         IssmDouble grade_g[NUMVERTICES] = {0.0};
    3329 
    3330         /* Get node coordinates and dof list: */
    3331         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3332         GradientIndexing(&vertexpidlist[0],control_index);
    3333 
    3334         /*Retrieve all inputs and parameters*/
    3335         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3336         this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    3337         this->parameters->FindParam(&responses,NULL,InversionCostFunctionsEnum);
    3338         Input* thickness_input            = inputs->GetInput(ThicknessEnum);                          _assert_(thickness_input);
    3339         Input* thicknessobs_input         = inputs->GetInput(InversionThicknessObsEnum);              _assert_(thicknessobs_input);
    3340         Input* weights_input              = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    3341         Input* vx_input                   = inputs->GetInput(VxEnum);                                 _assert_(vx_input);
    3342         Input* vy_input                   = inputs->GetInput(VyEnum);                                 _assert_(vy_input);
    3343         Input* surface_mass_balance_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum);         _assert_(surface_mass_balance_input);
    3344         Input* basal_melting_input        = inputs->GetInput(BasalforcingsMeltingRateEnum);           _assert_(basal_melting_input);
    3345         Input* dhdt_input                 = inputs->GetInput(BalancethicknessThickeningRateEnum);     _assert_(dhdt_input);
    3346 
    3347         /* Start  looping on the number of gaussian points: */
    3348         gauss=new GaussTria(2);
    3349         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3350 
    3351                 gauss->GaussPoint(ig);
    3352 
    3353                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3354                 GetNodalFunctions(basis, gauss);
    3355                 GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
    3356 
    3357                 thickness_input->GetInputValue(&thickness, gauss);
    3358                 thickness_input->GetInputDerivativeValue(&dH[0],&xyz_list[0][0],gauss);
    3359                 thicknessobs_input->GetInputValue(&thicknessobs, gauss);
    3360 
    3361                 /*Loop over all requested responses*/
    3362                 for(resp=0;resp<num_responses;resp++){
    3363 
    3364                         weights_input->GetInputValue(&weight,gauss,responses[resp]);
    3365 
    3366                         switch(responses[resp]){
    3367 
    3368                                 case ThicknessAbsMisfitEnum:
    3369                                         for(i=0;i<NUMVERTICES;i++) grade_g[i]+= (thicknessobs-thickness)*weight*Jdet*gauss->weight*basis[i];
    3370                                         break;
    3371                                 case ThicknessAbsGradientEnum:
    3372                                         for(i=0;i<NUMVERTICES;i++) grade_g[i]+= - weight*dH[0]*dbasis[0][i]*Jdet*gauss->weight;
    3373                                         for(i=0;i<NUMVERTICES;i++) grade_g[i]+= - weight*dH[1]*dbasis[1][i]*Jdet*gauss->weight;
    3374                                         break;
    3375                                 case ThicknessAlongGradientEnum:
    3376                                         vx_input->GetInputValue(&vx,gauss);
    3377                                         vy_input->GetInputValue(&vy,gauss);
    3378                                         vel = sqrt(vx*vx+vy*vy);
    3379                                         vx  = vx/(vel+1.e-9);
    3380                                         vy  = vy/(vel+1.e-9);
    3381                                         for(i=0;i<NUMVERTICES;i++) grade_g[i]+= - weight*(dH[0]*vx+dH[1]*vy)*(dbasis[0][i]*vx+dbasis[1][i]*vy)*Jdet*gauss->weight;
    3382                                         break;
    3383                                 case ThicknessAcrossGradientEnum:
    3384                                         vx_input->GetInputValue(&vx,gauss);
    3385                                         vy_input->GetInputValue(&vy,gauss);
    3386                                         vel = sqrt(vx*vx+vy*vy);
    3387                                         vx  = vx/(vel+1.e-9);
    3388                                         vy  = vy/(vel+1.e-9);
    3389                                         for(i=0;i<NUMVERTICES;i++) grade_g[i]+= - weight*(dH[0]*(-vy)+dH[1]*vx)*(dbasis[0][i]*(-vy)+dbasis[1][i]*vx)*Jdet*gauss->weight;
    3390                                         break;
    3391                                 case BalancethicknessMisfitEnum:
    3392                                         surface_mass_balance_input->GetInputValue(&surface_mass_balance,gauss);
    3393                                         basal_melting_input->GetInputValue(&basal_melting,gauss);
    3394                                         dhdt_input->GetInputValue(&dhdt,gauss);
    3395                                         vx_input->GetInputValue(&vx,gauss);
    3396                                         vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
    3397                                         vy_input->GetInputValue(&vy,gauss);
    3398                                         vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
    3399                                         for(i=0;i<NUMVERTICES;i++){
    3400                                                 grade_g[i]+= - weight*Jdet*gauss->weight*(
    3401                                                                         (vx*dH[0]+vy*dH[1] + thickness*(dvx[0]+dvy[1]))*(vx*dbasis[0][i]+ vy*dbasis[1][i] + basis[i]*(dvx[0]+dvy[1]))
    3402                                                                         -(surface_mass_balance-basal_melting-dhdt)*(vx*dbasis[0][i]+ vy*dbasis[1][i] + basis[i]*(dvx[0]+dvy[1]))
    3403                                                                         );
    3404                                         }
    3405                                         break;
    3406                                 default:
    3407                                         _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
    3408                         }
    3409                 }
    3410         }
    3411 
    3412         gradient->SetValues(NUMVERTICES,vertexpidlist,grade_g,ADD_VAL);
    3413 
    3414         /*Clean up and return*/
    3415         delete gauss;
    3416         xDelete<int>(responses);
    3417 }
    3418 /*}}}*/
    34192775void       Tria::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data){/*{{{*/
    34202776

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

    r18060 r18062  
    136136                #endif
    137137
    138                 void       Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index);
    139                 void       GradjBGradient(Vector<IssmDouble>* gradient,int control_index);
    140                 void       GradjDGradient(Vector<IssmDouble>* gradient,int control_index);
    141                 void       GradjBSSA(Vector<IssmDouble>* gradient,int control_index);
    142                 void       GradjDSSA(Vector<IssmDouble>* gradient,int control_index);
    143                 void       GradjDragSSA(Vector<IssmDouble>* gradient,int control_index);
    144                 void       GradjDragFS(Vector<IssmDouble>* gradient,int control_index);
    145                 void       GradjDragGradient(Vector<IssmDouble>* gradient,int control_index);
    146                 void       GradjDhDtBalancedthickness(Vector<IssmDouble>* gradient,int control_index);
    147                 void       GradjVxBalancedthickness(Vector<IssmDouble>* gradient,int control_index);
    148                 void       GradjVyBalancedthickness(Vector<IssmDouble>* gradient,int control_index);
    149                 void       GradjThicknessBalancethicknessSoft(Vector<IssmDouble>* gradient,int control_index);
    150                 void       GradjAdotBulancethickness2(Vector<IssmDouble>* gradient,int control_index);
    151138                void       GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data);
    152139                void       SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index);
Note: See TracChangeset for help on using the changeset viewer.