Changeset 6131


Ignore:
Timestamp:
10/04/10 08:00:41 (14 years ago)
Author:
Mathieu Morlighem
Message:

Added CM on vx and vy in CM balanced thickness

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

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/UpdateElementsBalancedthickness.cpp

    r5358 r6131  
    3939                IoModelFetchData(&iomodel->thickness_obs,NULL,NULL,iomodel_handle,"thickness_obs");
    4040                IoModelFetchData(&iomodel->weights,NULL,NULL,iomodel_handle,"weights");
    41                 if(iomodel->control_type==DhDtEnum){
     41                if(iomodel->control_type==DhDtEnum || iomodel->control_type==VxEnum || iomodel->control_type==VyEnum){
    4242                        IoModelFetchData(&iomodel->control_parameter,NULL,NULL,iomodel_handle,EnumToModelField(iomodel->control_type)); //copy the control parameter in iomodel
    4343                }

  • issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp

    r5631 r6131  
    2626
    2727                /*What control type?*/
    28                 if (iomodel->control_type!=DragCoefficientEnum && iomodel->control_type!=RheologyBbarEnum &&  iomodel->control_type!=DhDtEnum){
     28                if (iomodel->control_type!=DragCoefficientEnum &&
     29                                        iomodel->control_type!=RheologyBbarEnum &&
     30                                        iomodel->control_type!=DhDtEnum &&
     31                                        iomodel->control_type!=VxEnum &&
     32                                        iomodel->control_type!=VyEnum
     33                                        ){
    2934                        ISSMERROR("control_type %s not supported yet!",EnumToString(iomodel->control_type));
    3035                }

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

    r6130 r6131  
    673673                }
    674674                else if (control_type==RheologyBbarEnum){
    675                         Input* Bbar_input=matice->inputs->GetInput(RheologyBbarEnum); ISSMASSERT(Bbar_input);
    676                         Bbar_input->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], gauss);
    677                         //Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss->weight;
     675                        //nothing
    678676                }
    679677                else if (control_type==DhDtEnum){
    680                         Input* dhdt_input=inputs->GetInput(DhDtEnum);                 ISSMASSERT(dhdt_input);
    681                         dhdt_input->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], gauss);
    682                         //Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss->weight;
     678                        //nothing
     679                }
     680                else if (control_type==VxEnum){
     681                        //nothing
     682                }
     683                else if (control_type==VyEnum){
     684                        //nothing
    683685                }
    684686                else{
    685                         ISSMERROR("%s%i","unsupported control type: ",control_type);
     687                        ISSMERROR("unsupported control type: %s",EnumToString(control_type));
    686688                }
    687689        }
    … …  
    882884                        GradjDhDt(gradient);
    883885                        break;
     886                case VxEnum:
     887                        GradjVx(gradient);
     888                        break;
     889                case VyEnum:
     890                        GradjVy(gradient);
     891                        break;
    884892                default:
    885893                        ISSMERROR("%s%i","control type not supported yet: ",control_type);
    … …  
    964972
    965973        int i;
    966 
    967         /* node data: */
     974        int     ig;
     975
    968976        double       xyz_list[NUMVERTICES][3];
    969977        int          doflist1[NUMVERTICES];
    970978        double       dh1dh3[NDOF2][NUMVERTICES];
    971 
    972         /* gaussian points: */
    973         int     ig;
    974979        GaussTria *gauss=NULL;
    975 
    976         /* parameters: */
    977980        double  dk[NDOF2];
    978981        double  vx,vy;
    … …  
    983986        double  drag;
    984987        Friction* friction=NULL;
    985 
    986         /*element vector at the gaussian points: */
    987988        double  grade_g[NUMVERTICES]={0.0};
    988989        double  grade_g_gaussian[NUMVERTICES];
    989 
    990         /* Jacobian: */
    991990        double Jdet;
    992 
    993         /*nodal functions: */
    994991        double l1l2l3[3];
    995 
    996         /* strain rate: */
    997992        double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
    998 
    999         /*parameters: */
    1000993        double  cm_noisedmp;
    1001994        int analysis_type;
    … …  
    11001093}
    11011094/*}}}*/
     1095/*FUNCTION Tria::GradjVx{{{1*/
     1096void  Tria::GradjVx(Vec gradient){
     1097
     1098        /*Intermediaries*/
     1099        int        i,ig;
     1100        int        doflist1[NUMVERTICES];
     1101        double     l1l2l3[3];
     1102        double     thickness,Jdet;
     1103        double     Dlambda[2];
     1104        double     xyz_list[NUMVERTICES][3];
     1105        GaussTria *gauss                = NULL;
     1106        double     grade_g[NUMVERTICES] = {0.0};
     1107
     1108        /* Get node coordinates and dof list: */
     1109        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     1110        GetDofList1(&doflist1[0]);
     1111
     1112        /*Retrieve all inputs we will be needing: */
     1113        Input* adjoint_input=inputs->GetInput(AdjointEnum);     ISSMASSERT(adjoint_input);
     1114        Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input);
     1115
     1116        /* Start  looping on the number of gaussian points: */
     1117        gauss=new GaussTria(2);
     1118        for (ig=gauss->begin();ig<gauss->end();ig++){
     1119
     1120                gauss->GaussPoint(ig);
     1121
     1122                adjoint_input->GetParameterDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss);
     1123                thickness_input->GetParameterValue(&thickness, gauss);
     1124
     1125                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     1126                GetNodalFunctions(l1l2l3, gauss);
     1127
     1128                for(i=0;i<NUMVERTICES;i++) grade_g[i]=thickness*Dlambda[0]*Jdet*gauss->weight*l1l2l3[i];
     1129        }
     1130
     1131        VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);
     1132
     1133        /*Clean up and return*/
     1134        delete gauss;
     1135}
     1136/*}}}*/
     1137/*FUNCTION Tria::GradjVy{{{1*/
     1138void  Tria::GradjVy(Vec gradient){
     1139
     1140        /*Intermediaries*/
     1141        int        i,ig;
     1142        int        doflist1[NUMVERTICES];
     1143        double     thickness,Jdet;
     1144        double     l1l2l3[3];
     1145        double     Dlambda[2];
     1146        double     xyz_list[NUMVERTICES][3];
     1147        GaussTria *gauss                = NULL;
     1148        double     grade_g[NUMVERTICES] = {0.0};
     1149
     1150        /* Get node coordinates and dof list: */
     1151        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     1152        GetDofList1(&doflist1[0]);
     1153
     1154        /*Retrieve all inputs we will be needing: */
     1155        Input* adjoint_input=inputs->GetInput(AdjointEnum);     ISSMASSERT(adjoint_input);
     1156        Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input);
     1157
     1158        /* Start  looping on the number of gaussian points: */
     1159        gauss=new GaussTria(2);
     1160        for (ig=gauss->begin();ig<gauss->end();ig++){
     1161
     1162                gauss->GaussPoint(ig);
     1163
     1164                adjoint_input->GetParameterDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss);
     1165                thickness_input->GetParameterValue(&thickness, gauss);
     1166
     1167                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     1168                GetNodalFunctions(l1l2l3, gauss);
     1169
     1170                for(i=0;i<NUMVERTICES;i++) grade_g[i]=thickness*Dlambda[1]*Jdet*gauss->weight*l1l2l3[i];
     1171        }
     1172
     1173        VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);
     1174
     1175        /*Clean up and return*/
     1176        delete gauss;
     1177}
     1178/*}}}*/
    11021179/*FUNCTION Tria::InputControlUpdate{{{1*/
    11031180void  Tria::InputControlUpdate(double scalar,bool save_parameter){
    … …  
    11141191
    11151192        /*Rheology*/
    1116         if(control_type==RheologyBbarEnum){
    1117 
    1118                 /*First, get revert to previous parameter value (kept in ControlParameter input)*/
    1119                 matice->inputs->DuplicateInput(ControlParameterEnum,RheologyBbarEnum);
    1120 
    1121                 /*Now, update using the gradient new = old + scalar * gradient*/
    1122                 //matice->inputs->AXPY(RheologyBbarEnum,scalar,GradientEnum);
    1123                 // For now: Gradient is in element (TO BE CHANGED) and parameter in matice
    1124                 Input* input_B   =(Input*)matice->inputs->GetInput(RheologyBbarEnum); ISSMASSERT(input_B);
    1125                 Input* input_Grad=(Input*)this->inputs->GetInput(GradientEnum); ISSMASSERT(input_Grad);
    1126                 input_B->AXPY(input_Grad,scalar);
    1127 
    1128                 /*Constrain constrain input*/
    1129                 input=(Input*)matice->inputs->GetInput(RheologyBbarEnum); ISSMASSERT(input);
    1130                 input->Constrain(cm_min,cm_max);
    1131 
    1132                 /*Finally: save input if requested*/
    1133                 if (save_parameter) matice->inputs->DuplicateInput(RheologyBbarEnum,ControlParameterEnum);
    1134 
    1135         }
    1136         else if(control_type==DragCoefficientEnum){
    1137 
    1138                 /*First, get revert to previous parameter value (kept in ControlParameter input)*/
    1139                 this->inputs->DuplicateInput(ControlParameterEnum,DragCoefficientEnum);
    1140 
    1141                 /*Now, update using the gradient new = old + scalar * gradient*/
    1142                 this->inputs->AXPY(DragCoefficientEnum,scalar,GradientEnum);
    1143 
    1144                 /*Constrain input*/
    1145                 input=(Input*)this->inputs->GetInput(DragCoefficientEnum); ISSMASSERT(input);
    1146                 input->Constrain(cm_min,cm_max);
    1147 
    1148                 /*Finally: save input if requested*/
    1149                 if (save_parameter) inputs->DuplicateInput(DragCoefficientEnum,ControlParameterEnum);
    1150 
    1151         }
    1152         else if(control_type==DhDtEnum){
    1153 
    1154                 /*First, get revert to previous parameter value (kept in ControlParameter input)*/
    1155                 this->inputs->DuplicateInput(ControlParameterEnum,DhDtEnum);
    1156 
    1157                 /*Now, update using the gradient new = old + scalar * gradient*/
    1158                 this->inputs->AXPY(DhDtEnum,scalar,GradientEnum);
    1159 
    1160                 /*Constrain input*/
    1161                 input=(Input*)this->inputs->GetInput(DhDtEnum); ISSMASSERT(input);
    1162                 input->Constrain(cm_min,cm_max);
    1163 
    1164                 /*Finally: save input if requested*/
    1165                 if (save_parameter) inputs->DuplicateInput(DhDtEnum,ControlParameterEnum);
    1166 
    1167         }
    1168         else{
    1169                 ISSMERROR("control type %s not implemented yet",EnumToString(control_type));
     1193        switch(control_type){
     1194
     1195                case RheologyBbarEnum:{
     1196
     1197                        /*First, get revert to previous parameter value (kept in ControlParameter input)*/
     1198                        matice->inputs->DuplicateInput(ControlParameterEnum,RheologyBbarEnum);
     1199
     1200                        /*Now, update using the gradient new = old + scalar * gradient*/
     1201                        //matice->inputs->AXPY(RheologyBbarEnum,scalar,GradientEnum);
     1202                        // For now: Gradient is in element (TO BE CHANGED) and parameter in matice
     1203                        Input* input_B   =(Input*)matice->inputs->GetInput(RheologyBbarEnum); ISSMASSERT(input_B);
     1204                        Input* input_Grad=(Input*)this->inputs->GetInput(GradientEnum); ISSMASSERT(input_Grad);
     1205                        input_B->AXPY(input_Grad,scalar);
     1206
     1207                        /*Constrain constrain input*/
     1208                        input=(Input*)matice->inputs->GetInput(RheologyBbarEnum); ISSMASSERT(input);
     1209                        input->Constrain(cm_min,cm_max);
     1210
     1211                        /*Finally: save input if requested*/
     1212                        if (save_parameter) matice->inputs->DuplicateInput(RheologyBbarEnum,ControlParameterEnum);
     1213                                                                         }
     1214                        break;
     1215
     1216                case DragCoefficientEnum:{
     1217                        this->inputs->DuplicateInput(ControlParameterEnum,DragCoefficientEnum);
     1218                        this->inputs->AXPY(DragCoefficientEnum,scalar,GradientEnum);
     1219                        input=(Input*)this->inputs->GetInput(DragCoefficientEnum); ISSMASSERT(input);
     1220                        input->Constrain(cm_min,cm_max);
     1221                        if (save_parameter) inputs->DuplicateInput(DragCoefficientEnum,ControlParameterEnum);
     1222                                                                                 }
     1223                        break;
     1224                case DhDtEnum:{
     1225                        this->inputs->DuplicateInput(ControlParameterEnum,DhDtEnum);
     1226                        this->inputs->AXPY(DhDtEnum,scalar,GradientEnum);
     1227                        input=(Input*)this->inputs->GetInput(DhDtEnum); ISSMASSERT(input);
     1228                        input->Constrain(cm_min,cm_max);
     1229                        if (save_parameter) inputs->DuplicateInput(DhDtEnum,ControlParameterEnum);
     1230                                                  }
     1231                        break;
     1232                case VxEnum:{
     1233                        this->inputs->DuplicateInput(ControlParameterEnum,VxEnum);
     1234                        this->inputs->AXPY(VxEnum,scalar,GradientEnum);
     1235                        input=(Input*)this->inputs->GetInput(VxEnum); ISSMASSERT(input);
     1236                        input->Constrain(cm_min,cm_max);
     1237                        if (save_parameter) inputs->DuplicateInput(VxEnum,ControlParameterEnum);
     1238                                                  }
     1239                        break;
     1240                case VyEnum:{
     1241                        this->inputs->DuplicateInput(ControlParameterEnum,VyEnum);
     1242                        this->inputs->AXPY(VyEnum,scalar,GradientEnum);
     1243                        input=(Input*)this->inputs->GetInput(VyEnum); ISSMASSERT(input);
     1244                        input->Constrain(cm_min,cm_max);
     1245                        if (save_parameter) inputs->DuplicateInput(VyEnum,ControlParameterEnum);
     1246                                                  }
     1247                        break;
     1248                default:
     1249                        ISSMERROR("control type %s not implemented yet",EnumToString(control_type));
    11701250        }
    11711251
    … …  
    38943974                 *      2  [       obs            obs   ]
    38953975                 *
    3896                  *        dJ             2
     3976                 *        dJ
    38973977                 * DU = - -- = (u   - u )
    38983978                 *        du     obs

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

    r6130 r6131  
    8686                void   GradjDrag(Vec gradient);
    8787                void   GradjDhDt(Vec gradient);
     88                void   GradjVx(Vec gradient);
     89                void   GradjVy(Vec gradient);
    8890                void   InputControlUpdate(double scalar,bool save_parameter);
    8991                void   InputArtificialNoise(int enum_type,double min, double max);

  • issm/trunk/src/c/solutions/adjointbalancedthickness_core.cpp

    r5414 r6131  
    3838       
    3939        /*Save results*/
    40         if(solution_type==AdjointSolutionEnum && !control_analysis){
     40        //if(solution_type==AdjointSolutionEnum && !control_analysis){
    4141                if(verbose)_printf_("saving results:\n");
    4242                InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointEnum);
    43         }
     43        //}
    4444}
Note: See TracChangeset for help on using the changeset viewer.