Changeset 22739


Ignore:
Timestamp:
05/03/18 13:31:36 (7 years ago)
Author:
erobo
Message:

CHG: allow for transient ControlInputs

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

Legend:

Unmodified
Added
Removed

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

    r22625 r22739  
    13871387        /*Get number of controls*/
    13881388        int num_controls;
     1389        bool isautodiff;
    13891390        parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     1391        parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
     1392
    13901393
    13911394        /*Get number of vertices*/
    13921395        int numvertices = this->GetNumberOfVertices();
    1393 
     1396        if(isautodiff){
     1397                int* N;
     1398                int* M;
     1399                int start = 0;
     1400                parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
     1401                parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
     1402                if(control_index>0) for(int n=0;n<control_index-1;n++) start+=N[n]*M[n];
     1403
     1404                for(int n=0;n<N[control_index];n++){
     1405                        for(int i=0;i<numvertices;i++){
     1406                                indexing[i+n*numvertices]=this->vertices[i]->Sid() + start + n*M[control_index];
     1407                        }
     1408                }
     1409        }
     1410        else{
     1411                int M;
     1412                parameters->FindParam(&M,ControlInputSizeMEnum);
    13941413        /*get gradient indices*/
    1395         if(onsid){
    1396                 for(int i=0;i<numvertices;i++){
    1397                         indexing[i]=num_controls*this->vertices[i]->Sid() + control_index;
    1398                 }
    1399         }
    1400         else{
    1401                 for(int i=0;i<numvertices;i++){
    1402                         indexing[i]=num_controls*this->vertices[i]->Pid() + control_index;
    1403                 }
    1404         }
    1405 
     1414                if(onsid){
     1415                        for(int i=0;i<numvertices;i++){
     1416                                indexing[i]=this->vertices[i]->Sid() + (control_index)*M;
     1417                        }
     1418                }
     1419                else{
     1420                        for(int i=0;i<numvertices;i++){
     1421                                indexing[i]=this->vertices[i]->Pid() + (control_index)*M;
     1422                        }
     1423                }
     1424        }
    14061425}
    14071426/*}}}*/
    … …  
    16191638
    16201639        /*For now we only support nodal vectors*/
    1621         if(M!=iomodel->numberofvertices) _error_("not supported");
    1622         if(N!=1) _error_("not supported");
     1640        //if(M!=iomodel->numberofvertices) _error_("not supported");
     1641        //if(N!=1) _error_("not supported");
    16231642
    16241643        /*Recover vertices ids needed to initialize inputs*/
    … …  
    16371656                this->AddControlInput(input_enum,values,min_vector,max_vector,P1Enum,id);
    16381657        }
     1658
     1659        else if(M==iomodel->numberofvertices+1){
     1660            /*create transient input: */
     1661            IssmDouble* times = xNew<IssmDouble>(N);
     1662            for(int t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
     1663                                /*Create the three transient inputs for the control input*/
     1664            TransientInput* values_input=new TransientInput(ControlInputValuesEnum,times,N);
     1665                                TransientInput* mins_input = new TransientInput(ControlInputMinsEnum,times,N);
     1666                                TransientInput* maxs_input = new TransientInput(ControlInputMaxsEnum,times,N);
     1667                           TransientInput* grad_input = new TransientInput(ControlInputGradEnum);
     1668                                for(int t=0;t<N;t++){
     1669                for(int i=0;i<numvertices;i++){
     1670                                                values[i]=vector[N*(vertexids[i]-1)+t];
     1671                                                values_min[i] = min_vector[N*(vertexids[i]-1)+t];
     1672                                                values_max[i] = max_vector[N*(vertexids[i]-1)+t];
     1673                                         }
     1674                                        switch(this->ObjectEnum()){
     1675                    case TriaEnum:
     1676                                                                        values_input->AddTimeInput(new TriaInput(ControlInputValuesEnum,values,P1Enum));
     1677                                                                        mins_input->AddTimeInput(new TriaInput(ControlInputMinsEnum,values_min,P1Enum));
     1678                                                                        maxs_input->AddTimeInput(new TriaInput(ControlInputMaxsEnum,values_max,P1Enum));
     1679                                                                break;
     1680                    case PentaEnum:
     1681                                                                printf("-------------- file: Element.cpp line: %i\n",__LINE__);
     1682                                                                        values_input->AddTimeInput(new PentaInput(ControlInputValuesEnum,values,P1Enum));
     1683                                                                        mins_input->AddTimeInput(new PentaInput(ControlInputMinsEnum,values_min,P1Enum));
     1684                                                                        maxs_input->AddTimeInput(new PentaInput(ControlInputMaxsEnum,values_max,P1Enum));
     1685                                                                        break;
     1686                    case TetraEnum:
     1687                                                                        values_input->AddTimeInput(new TetraInput(ControlInputValuesEnum,values,P1Enum));
     1688                                                                        mins_input->AddTimeInput(new TetraInput(ControlInputMinsEnum,values_min,P1Enum));
     1689                                                                        maxs_input->AddTimeInput(new TetraInput(ControlInputMaxsEnum,values_max,P1Enum));
     1690                                                                        break;
     1691                    default: _error_("Not implemented yet");
     1692                }
     1693            }
     1694            this->inputs->AddInput(new ControlInput(input_enum,TransientInputEnum,values_input,mins_input,maxs_input,grad_input,P1Enum,id));
     1695            xDelete<IssmDouble>(times);
     1696        }
     1697                  else _error_("not currently supported type of M and N attempted");
    16391698
    16401699        /*clean up*/

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

    r22732 r22739  
    202202                virtual void       ComputeEsaStrainAndVorticity(void)=0;
    203203                virtual void       Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0;
     204                virtual void       ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N,int M)=0;
    204205                virtual void       ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index)=0;
    205206                virtual void       ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum)=0;
    … …  
    224225                virtual int        GetNumberOfNodes(int enum_type)=0;
    225226                virtual int        GetNumberOfVertices(void)=0;
    226                 virtual void       GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,bool onsid)=0;
     227                virtual void       GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,int offset, bool onsid)=0;
     228                virtual void       GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data, bool onsid)=0;
    227229                virtual void       GetVerticesCoordinatesBase(IssmDouble** xyz_list)=0;
    228230                virtual void       GetVerticesCoordinatesTop(IssmDouble** xyz_list)=0;
    … …  
    285287                virtual void       ResetHooks()=0;
    286288                virtual void       ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments){_error_("not implemented yet");};
     289                virtual void       SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int N,int M)=0;
    287290                virtual void       SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index)=0;
    288291                virtual void       SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters)=0;

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

    r22647 r22739  
    561561}
    562562/*}}}*/
     563void       Penta::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N,int M){/*{{{*/
     564
     565        int    vertexpidlist[NUMVERTICES];
     566        IssmDouble grad_list[NUMVERTICES];
     567        Input* grad_input=NULL;
     568        Input* input=NULL;
     569
     570        if(enum_type==MaterialsRheologyBbarEnum){
     571                input=(Input*)inputs->GetInput(MaterialsRheologyBEnum);
     572        }
     573        else if(enum_type==DamageDbarEnum){
     574                input=(Input*)inputs->GetInput(DamageDEnum);
     575        }
     576        else{
     577                input=inputs->GetInput(enum_type);
     578        }
     579        if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found");
     580        if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(enum_type) << " is not a ControlInput");
     581
     582        GradientIndexing(&vertexpidlist[0],control_index);
     583        for(int i=0;i<NUMVERTICES;i++) grad_list[i]=gradient[vertexpidlist[i]];
     584        grad_input=new PentaInput(GradientEnum,grad_list,P1Enum);
     585        ((ControlInput*)input)->SetGradient(grad_input);
     586
     587}
     588/*}}}*/
    563589void       Penta::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){/*{{{*/
    564590
    … …  
    11491175        ((ControlInput*)input)->GetVectorFromInputs(vector,&vertexidlist[0],data);
    11501176}
     1177/*}}}*/
     1178void       Penta::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data,int offset,bool onsid){/*{{{*/
     1179
     1180        int* idlist = NULL;
     1181        IssmDouble* values = NULL;
     1182        int* M;
     1183
     1184        /*Get out if this is not an element input*/
     1185        if(!IsInput(control_enum)) _error_("Enum "<<EnumToStringx(control_enum)<<" is not in IsInput");
     1186        Input* input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
     1187
     1188
     1189        /*Cast to Controlinput*/
     1190        if(input->ObjectEnum()!=ControlInputEnum) _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
     1191        ControlInput* controlinput = xDynamicCast<ControlInput*>(input);
     1192
     1193        if(strcmp(data,"value")==0){
     1194                input  = controlinput->values;
     1195        }
     1196        else if (strcmp(data,"lowerbound")==0){
     1197                input = controlinput->minvalues;
     1198        }
     1199        else if (strcmp(data,"upperbound")==0){
     1200                input = controlinput->maxvalues;
     1201        }
     1202        else if (strcmp(data,"gradient")==0){
     1203                input = controlinput->gradient;
     1204        }
     1205        else{
     1206                _error_("Data " << data << " not supported yet");
     1207        }
     1208        /*Check what input we are dealing with*/
     1209
     1210        switch(input->ObjectEnum()){
     1211                case PentaInputEnum:
     1212                                  {
     1213                                        PentaInput* pentainput = xDynamicCast<PentaInput*>(input);
     1214                                        if(pentainput->interpolation_type!=P1Enum) _error_("not supported yet");
     1215
     1216                                        /*Create list of indices and values for global vector*/
     1217                                        idlist = xNew<int>(NUMVERTICES);
     1218                                        values = xNew<IssmDouble>(NUMVERTICES);
     1219                                        GradientIndexing(&idlist[0],control_index,true);
     1220                                        for(int i=0;i<NUMVERTICES;i++){
     1221                                                values[i] = pentainput->values[i];
     1222                                                //if(this->Id()<=10)_printf_("index "<<idlist[i]<<"\n");
     1223                                        }
     1224                                        vector->SetValues(NUMVERTICES,idlist,values,INS_VAL);
     1225                                        break;
     1226                                  }
     1227
     1228                                case TransientInputEnum:
     1229                                  {
     1230                                        parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
     1231                                        TransientInput* transientinput = xDynamicCast<TransientInput*>(input);
     1232                                        int N = transientinput->numtimesteps;
     1233                                        idlist = xNew<int>(NUMVERTICES*N);
     1234                                        values = xNew<IssmDouble>(NUMVERTICES*N);
     1235                                        for(int t=0;t<transientinput->numtimesteps;t++) {
     1236                                                IssmDouble time = transientinput->GetTimeByOffset(t);
     1237                                                input = transientinput->GetTimeInput(time);
     1238                                                TriaInput* timeinput = xDynamicCast<TriaInput*>(input);
     1239                                                if(timeinput->interpolation_type!=P1Enum) _error_("not supported yet");
     1240                                                /*Create list of indices and values for global vector*/
     1241                                                for(int i=0;i<NUMVERTICES;i++){
     1242                                                        idlist[N*i+t] = offset + this->vertices[i]->Sid()+t*M[control_index];
     1243                                                        values[N*i+t] = timeinput->values[i];
     1244                                                }
     1245                                        }
     1246
     1247                                        vector->SetValues(NUMVERTICES*transientinput->numtimesteps,idlist,values,INS_VAL);
     1248                                        break;
     1249                                  }
     1250                                default: _error_("input "<<input->ObjectEnum()<<" not supported yet");
     1251
     1252                }
     1253                        xDelete<int>(idlist);
     1254                        xDelete<IssmDouble>(values);
     1255        }
    11511256/*}}}*/
    11521257void       Penta::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){/*{{{*/
    … …  
    22522357}
    22532358
     2359/*}}}*/
     2360void       Penta::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset, int N, int M){/*{{{*/
     2361
     2362        IssmDouble  values[NUMVERTICES];
     2363        int         vertexpidlist[NUMVERTICES],control_init;
     2364
     2365        /*Specific case for depth averaged quantities*/
     2366        control_init=control_enum;
     2367        if(control_enum==MaterialsRheologyBbarEnum){
     2368                control_enum=MaterialsRheologyBEnum;
     2369                if(!IsOnBase()) return;
     2370        }
     2371        if(control_enum==DamageDbarEnum){
     2372                control_enum=DamageDEnum;
     2373                if(!IsOnBase()) return;
     2374        }
     2375
     2376        /*Get out if this is not an element input*/
     2377        if(!IsInput(control_enum)) return;
     2378
     2379        /*Prepare index list*/
     2380        GradientIndexing(&vertexpidlist[0],control_index);
     2381
     2382        /*Get values on vertices*/
     2383        for(int i=0;i<NUMVERTICES;i++){
     2384                values[i]=vector[vertexpidlist[i]];
     2385        }
     2386        Input* new_input = new PentaInput(control_enum,values,P1Enum);
     2387        Input* input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
     2388        if(input->ObjectEnum()!=ControlInputEnum){
     2389                _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
     2390        }
     2391
     2392        ((ControlInput*)input)->SetInput(new_input);
     2393
     2394        if(control_init==MaterialsRheologyBbarEnum){
     2395                this->InputExtrude(control_enum,-1);
     2396        }
     2397        if(control_init==DamageDbarEnum){
     2398                this->InputExtrude(control_enum,-1);
     2399        }
     2400}
    22542401/*}}}*/
    22552402void       Penta::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){/*{{{*/

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

    r22732 r22739  
    5858                void           ComputeStressTensor();
    5959                void           Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
     60                void           ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N,int M);
    6061                void           ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index);
    6162                void           ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
    … …  
    8687                Penta*         GetUpperPenta(void);
    8788                void           GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,bool onsid);
     89                void           GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,int offset,bool onsid);
    8890                void           GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
    8991                void           GetVerticesCoordinatesTop(IssmDouble** pxyz_list);
    … …  
    148150                void           ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments);
    149151                void             SetClone(int* minranks);
     152                void           SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset, int N,int M);
    150153                void           SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index);
    151154                void           SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);

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

    r22732 r22739  
    5252                void        ComputeStressTensor(){_error_("not implemented yet");};
    5353                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters){_error_("not implemented yet");};
     54                void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N,int M){_error_("not implemented yet");};
    5455                void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");};
    5556                void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){_error_("not implemented yet");};
    … …  
    7273                int         GetNumberOfNodes(int enum_type){_error_("not implemented yet");};
    7374                int         GetNumberOfVertices(void);
     75                void        GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,int offset,bool onsid){_error_("not implemented yet");};
    7476                void        GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,bool onsid){_error_("not implemented yet");};
    7577                void        GetVerticesCoordinates(IssmDouble** pxyz_list);
    … …  
    135137                void        ResetFSBasalBoundaryCondition(void){_error_("not implemented yet");};
    136138                void        ResetHooks(){_error_("not implemented yet");};
     139                void        SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int N,int M){_error_("not implemented yet");};
    137140                void        SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){_error_("not implemented yet");};
    138141                void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters){_error_("not implemented yet");};

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

    r22732 r22739  
    5252                void        ComputeEsaStrainAndVorticity(){_error_("not implemented yet!");};
    5353                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
     54                void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N,int M){_error_("not implemented yet");};
    5455                void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");};
    5556                void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){_error_("not implemented yet");};
    … …  
    7879                int         GetNumberOfNodes(int enum_type){_error_("not implemented yet");};
    7980                int         GetNumberOfVertices(void);
     81                void        GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,int offset,bool onsid){_error_("not implemented yet");};
    8082                void        GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,bool onsid){_error_("not implemented yet");};
    8183                void        GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
    … …  
    143145                void        ResetHooks();
    144146                void        ReduceMatrices(ElementMatrix* Ke,ElementVector* pe);
     147                void        SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int N, int M){_error_("not implemented yet");};
    145148                void        SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){_error_("not implemented yet");};
    146149                void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);

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

    r22732 r22739  
    742742}
    743743/*}}}*/
     744void       Tria::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N, int M){/*{{{*/
     745
     746        int    idlist[NUMVERTICES];
     747        int     gradidlist[NUMVERTICES];
     748        IssmDouble grad_list[NUMVERTICES];
     749        Input* grad_input=NULL;
     750
     751        Input* input=inputs->GetInput(enum_type);
     752        if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found");
     753        if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(enum_type) << " is not a ControlInput");
     754
     755        GradientIndexing(&gradidlist[0],control_index,true);
     756
     757        for(int n=0;n<N;n++){
     758                for(int i=0;i<NUMVERTICES;i++){
     759                        idlist[i] = offset + this->vertices[i]->Sid()+n*M;
     760                        grad_list[i]=gradient[idlist[i]];
     761                }
     762
     763                ControlInput* controlinput = xDynamicCast<ControlInput*>(input);
     764                if(controlinput->layout_enum!=TransientInputEnum){
     765                        grad_input=new TriaInput(GradientEnum,grad_list,P1Enum);
     766                        controlinput->SetGradient(grad_input);
     767                }
     768                else{
     769                        grad_input = new TriaInput(GradientEnum,grad_list,P1Enum);
     770                        controlinput->SetGradient(grad_input,n);
     771                        controlinput->Configure(parameters);
     772                }
     773        }
     774
     775               
     776}/*}}}*/
    744777void       Tria::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){/*{{{*/
    745778
    … …  
    16761709
    16771710        ((ControlInput*)input)->GetVectorFromInputs(vector,&vertexidlist[0],data);
     1711}
     1712/*}}}*/
     1713void       Tria::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data,int offset, bool onsid){/*{{{*/
     1714
     1715        int* idlist = NULL;
     1716        IssmDouble* values = NULL;
     1717        int* M;
     1718
     1719        /*Get out if this is not an element input*/
     1720        if(!IsInput(control_enum)) _error_("Enum "<<EnumToStringx(control_enum)<<" is not in IsInput");
     1721        Input* input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
     1722
     1723        parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
     1724       
     1725        /*Cast to Controlinput*/
     1726        if(input->ObjectEnum()!=ControlInputEnum) _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
     1727        ControlInput* controlinput = xDynamicCast<ControlInput*>(input);
     1728       
     1729        if(strcmp(data,"value")==0){
     1730                input  = controlinput->values;
     1731        }
     1732        else if (strcmp(data,"lowerbound")==0){
     1733                input = controlinput->minvalues;
     1734        }
     1735        else if (strcmp(data,"upperbound")==0){
     1736                input = controlinput->maxvalues;
     1737        }
     1738        else if (strcmp(data,"gradient")==0){
     1739                input = controlinput->gradient;
     1740        }
     1741        else{
     1742                _error_("Data " << data << " not supported yet");
     1743        }
     1744        /*Check what input we are dealing with*/
     1745       
     1746        switch(input->ObjectEnum()){
     1747                case TriaInputEnum:
     1748                          {
     1749                                TriaInput* triainput = xDynamicCast<TriaInput*>(input);
     1750                                if(triainput->interpolation_type!=P1Enum) _error_("not supported yet");
     1751
     1752                                /*Create list of indices and values for global vector*/
     1753                                idlist = xNew<int>(NUMVERTICES);
     1754                                values = xNew<IssmDouble>(NUMVERTICES);
     1755                                GradientIndexing(&idlist[0],control_index,true);
     1756                                for(int i=0;i<NUMVERTICES;i++){
     1757                                        values[i] = triainput->values[i];
     1758                                }
     1759                                vector->SetValues(NUMVERTICES,idlist,values,INS_VAL);
     1760                                break;
     1761                          }
     1762                       
     1763                case TransientInputEnum:
     1764                                {
     1765                                        TransientInput* transientinput = xDynamicCast<TransientInput*>(input);
     1766                                        int N = transientinput->numtimesteps;
     1767                                        idlist = xNew<int>(NUMVERTICES*N);
     1768                                        values = xNew<IssmDouble>(NUMVERTICES*N);
     1769                                        for(int t=0;t<transientinput->numtimesteps;t++) {
     1770                                                IssmDouble time = transientinput->GetTimeByOffset(t);
     1771                                                input = transientinput->GetTimeInput(time);
     1772                                                TriaInput* timeinput = xDynamicCast<TriaInput*>(input);
     1773                                                if(timeinput->interpolation_type!=P1Enum) _error_("not supported yet");
     1774                                                /*Create list of indices and values for global vector*/
     1775                                                for(int i=0;i<NUMVERTICES;i++){
     1776                                                                idlist[N*i+t] = offset + this->vertices[i]->Sid()+t*M[control_index];
     1777                                                                values[N*i+t] = timeinput->values[i];
     1778                                                }
     1779                                        }
     1780                                        vector->SetValues(NUMVERTICES*transientinput->numtimesteps,idlist,values,INS_VAL);
     1781                                        break;
     1782                                }
     1783                default: _error_("input "<<input->ObjectEnum()<<" not supported yet");
     1784        }
     1785        xDelete<int>(idlist);
     1786        xDelete<IssmDouble>(values);
    16781787}
    16791788/*}}}*/
    … …  
    19522061
    19532062        /*Get parameters: */
    1954         iomodel->FindConstant(&yts,"md.constants.yts"); 
     2063        iomodel->FindConstant(&yts,"md.constants.yts");
    19552064        iomodel->FindConstant(&control_analysis,"md.inversion.iscontrol");
    19562065        iomodel->FindConstant(&ad_analysis, "md.autodiff.isautodiff");
    … …  
    19632072
    19642073        /*Recover vertices ids needed to initialize inputs*/
    1965         for(i=0;i<3;i++){ 
     2074        for(i=0;i<3;i++){
    19662075                tria_vertex_ids[i]=reCast<int>(iomodel->elements[3*index+i]); //ids for vertices are in the elements array from Matlab
    19672076        }
    … …  
    20632172        /*DatasetInputs*/
    20642173        if (control_analysis && iomodel->Data("md.inversion.cost_functions_coefficients")){
    2065        
     2174
    20662175                /*Generate cost functions associated with the iomodel*/
    2067                 char**  cost_functions                  = NULL;
    2068                 int*            cost_functions_enums = NULL;
    2069                 int             num_cost_functions;
     2176                char**   cost_functions       = NULL;
     2177                int*     cost_functions_enums = NULL;
     2178                int      num_cost_functions;
    20702179
    20712180                iomodel->FindConstant(&num_cost_functions,"md.inversion.num_cost_functions");
    … …  
    20912200        }
    20922201}
    2093 /*}}}*/
    20942202void       Tria::InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type){/*{{{*/
    20952203
    … …  
    31493257
    31503258        _error_("not implemented yet");
     3259}
     3260/*}}}*/
     3261void       Tria::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int N, int M){/*{{{*/
     3262
     3263        IssmDouble  values[NUMVERTICES];
     3264        int         idlist[NUMVERTICES],control_init;
     3265
     3266
     3267        /*Get Domain type*/
     3268        int domaintype;
     3269        parameters->FindParam(&domaintype,DomainTypeEnum);
     3270
     3271        /*Specific case for depth averaged quantities*/
     3272        control_init=control_enum;
     3273        if(domaintype==Domain2DverticalEnum){
     3274                if(control_enum==MaterialsRheologyBbarEnum){
     3275                        control_enum=MaterialsRheologyBEnum;
     3276                        if(!IsOnBase()) return;
     3277                }
     3278                if(control_enum==DamageDbarEnum){
     3279                        control_enum=DamageDEnum;
     3280                        if(!IsOnBase()) return;
     3281                }
     3282        }
     3283
     3284        /*Get out if this is not an element input*/
     3285        if(!IsInput(control_enum)) return;
     3286       
     3287        Input* input     = (Input*)this->inputs->GetInput(control_enum);   _assert_(input);
     3288        if(input->ObjectEnum()!=ControlInputEnum){
     3289                _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
     3290        }
     3291
     3292        ControlInput* controlinput = xDynamicCast<ControlInput*>(input);
     3293        input = controlinput->values;
     3294
     3295        /*Get values on vertices*/
     3296        for(int n=0;n<N;n++){
     3297                for(int i=0;i<NUMVERTICES;i++){
     3298                        idlist[i]=offset + this->vertices[i]->Sid()+n*M;
     3299                        values[i]=vector[idlist[i]];
     3300                        //_printf_("index: "<<idlist[i]<<" -- value: "<<values[i]<<"\n");
     3301                }
     3302                if(input->ObjectEnum()==TriaInputEnum){
     3303                        Input* new_input = new TriaInput(control_enum,values,P1Enum);
     3304                        controlinput->SetInput(new_input);
     3305                }
     3306                else if(input->ObjectEnum()==TransientInputEnum){
     3307                        Input* new_input = new TriaInput(control_enum,values,P1Enum);
     3308                        controlinput->SetInput(new_input,n);
     3309                        controlinput->Configure(parameters);
     3310                }
     3311                else _error_("Type not supported");
     3312        }
    31513313}
    31523314/*}}}*/

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

    r22732 r22739  
    6363                void        ComputeSurfaceNormalVelocity();
    6464                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
     65                void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N, int M);
    6566                void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index);
    6667                void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
    … …  
    8586                int         GetNumberOfNodes(int enum_type);
    8687                int         GetNumberOfVertices(void);
     88                void        GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,int offset,bool onsid);
    8789                void        GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,bool onsid);
    8890                void        GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
    … …  
    118120                void        ResetHooks();
    119121                void        ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments);
     122                void        SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int N,int M);
    120123                void        SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index);
    121124                void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);

  • issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp

    r22515 r22739  
    2626        control_id=id;
    2727        enum_type=in_enum_type;
     28        layout_enum = input_layout_enum;
    2829
    2930        _assert_(interp==P1Enum);
    … …  
    4849}
    4950/*}}}*/
     51ControlInput::ControlInput(int in_enum_type,int input_layout_enum,Input* input_pvalues,Input* input_pmin,Input* input_pmax,Input* input_pgrad,int interp,int id){/*{{{*/
     52
     53        this->control_id=id;
     54        this->enum_type=in_enum_type;
     55        this->layout_enum = input_layout_enum;
     56
     57        _assert_(interp==P1Enum);
     58        if(input_layout_enum!=TransientInputEnum) _error_("Wrong type of layout_enum, needs to be a TransientInputEnum");
     59
     60                        this->values = input_pvalues;
     61                        this->savedvalues= NULL;
     62                        this->minvalues  = input_pmin;
     63                        this->maxvalues  = input_pmax;
     64                        this->gradient =  input_pgrad;
     65        }
     66/*}}}*/
    5067ControlInput::~ControlInput(){/*{{{*/
    51         delete values;
    52         delete savedvalues;
    53         delete minvalues;
    54         delete maxvalues;
    55         delete gradient;
     68        if(values)      delete values;
     69        if(savedvalues) delete savedvalues;
     70        if(minvalues)   delete minvalues;
     71        if(maxvalues)   delete maxvalues;
     72        if(gradient)    delete gradient;
    5673}
    5774/*}}}*/
    … …  
    6582        output->enum_type=this->enum_type;
    6683        output->control_id=this->control_id;
     84        output->layout_enum = this->control_id;
    6785
    6886        if(values)      output->values      = xDynamicCast<Input*>(this->values->copy());
    … …  
    7997        _printf_("ControlInput:\n");
    8098        _printf_(setw(15)<<"   ControlInput "<<setw(25)<<left<<EnumToStringx(this->enum_type)<<"\n");
     99        _printf_(setw(15)<<"   ControlInput "<<setw(25)<<left<<EnumToStringx(this->layout_enum)<<"\n");
    81100        _printf_("---values: \n");     if (values)      values->Echo();
    82101        _printf_("---savedvalues: \n");if (savedvalues) savedvalues->Echo();
    … …  
    98117        MARSHALLING(enum_type);
    99118        MARSHALLING(control_id);
     119        MARSHALLING(layout_enum);
    100120
    101121        if (marshall_direction == MARSHALLING_BACKWARD){
    … …  
    115135                                gradient   =new PentaInput();
    116136                                break;
     137                        case TransientInputEnum:
     138                                values    =new TransientInput();
     139                                savedvalues=new TransientInput();
     140                                minvalues  =new TransientInput();
     141                                maxvalues  =new TransientInput();
     142                                gradient   =new TransientInput();
     143                                break;
    117144                        default:
    118145                                _error_("Input of Enum " << EnumToStringx(enum_type) << " not supported yet by ControlInput");
    … …  
    146173}/*}}}*/
    147174void ControlInput::Configure(Parameters* parameters){/*{{{*/
     175        if(this->values->ObjectEnum()==TransientInputEnum){
     176                this->values->Configure(parameters);
     177                this->minvalues->Configure(parameters);
     178                this->maxvalues->Configure(parameters);
     179                this->gradient->Configure(parameters);
     180        }
    148181        /*do nothing: */
    149182}
    … …  
    194227}/*}}}*/
    195228void ControlInput::GetVectorFromInputs(Vector<IssmDouble>* vector,int* doflist,const char* data){/*{{{*/
    196          if(strcmp(data,"value")==0){
     229        if(strcmp(data,"value")==0){
    197230                 _assert_(values);
    198231                 values->GetVectorFromInputs(vector,doflist);
    … …  
    226259        this->savedvalues=xDynamicCast<Input*>(this->values->copy());
    227260}/*}}}*/
    228 void ControlInput::SetGradient(Input* gradient_in){/*{{{*/
     261void ControlInput::SetGradient(Input* gradient_in,int timestep){/*{{{*/
     262if(this->values->ObjectEnum()!=TransientInputEnum)_error_("you are in the wrong place, go home");
    229263
    230264        /*Get enum for current gradient*/
    … …  
    243277        }
    244278
     279        TransientInput* transient_input = xDynamicCast<TransientInput*>(this->gradient);
     280        TransientInput* values_input = xDynamicCast<TransientInput*>(this->values);
     281        if(values_input->numtimesteps==transient_input->numtimesteps){
     282                TransientInput* new_trans_input = new TransientInput(ControlInputGradEnum);
     283                IssmDouble time = transient_input->GetTimeByOffset(timestep);
     284                for(int i=0;i<transient_input->numtimesteps;i++){
     285                        if(transient_input->timesteps[i]==time) new_trans_input->AddTimeInput(xDynamicCast<TriaInput*>(gradient_in),time);
     286                        else {
     287                                Input* input = transient_input->GetTimeInput(transient_input->timesteps[i]);
     288                                new_trans_input->AddTimeInput(xDynamicCast<TriaInput*>(input),transient_input->timesteps[i]);
     289                        }
     290                }
     291                this->gradient=new_trans_input;
     292        }
     293        else{
     294                IssmDouble time = values_input->GetTimeByOffset(timestep);
     295                transient_input->AddTimeInput(gradient_in,time);
     296        }
     297}/*}}}*/
     298void ControlInput::SetGradient(Input* gradient_in){/*{{{*/
     299
     300        /*Get enum for current gradient*/
     301        switch(this->control_id){
     302                case 1:
     303                        gradient_in->ChangeEnum(Gradient1Enum);
     304                        break;
     305                case 2:
     306                        gradient_in->ChangeEnum(Gradient2Enum);
     307                        break;
     308                case 3:
     309                        gradient_in->ChangeEnum(Gradient3Enum);
     310                        break;
     311                default:
     312                        _error_("more than 3 controls not implemented yet (Gradient " << this->control_id << " was requested). EnumDefinitions.h needs to be updated.");
     313        }
     314
    245315        /*Delete old gradient and assign new gradient*/
    246316        if(gradient) delete gradient;
    … …  
    250320void ControlInput::SetInput(Input* in_input){/*{{{*/
    251321
     322        if(layout_enum==TransientInputEnum)_error_("need two arguments in SetInput for TransientInput Controls");
    252323        delete values; this->values=in_input;
    253324        this->SaveValue(); //because this is what SpawnResult saves FIXME
     325
     326}/*}}}*/
     327void ControlInput::SetInput(Input* in_input,int timeoffset){/*{{{*/
     328        Input* input = this->values;
     329        if(input->ObjectEnum()!=TransientInputEnum)_error_("cannot have timeoffset argument if not TransientInput Control");
     330        TransientInput* transient_input = xDynamicCast<TransientInput*>(input);
     331        IssmDouble time = transient_input->GetTimeByOffset(timeoffset);
     332        TransientInput* new_trans_input = new TransientInput(ControlInputValuesEnum);
     333        for(int i=0;i<transient_input->numtimesteps;i++){
     334                if(transient_input->timesteps[i]==time) new_trans_input->AddTimeInput(xDynamicCast<TriaInput*>(in_input),time);
     335                else {
     336                        input = transient_input->GetTimeInput(transient_input->timesteps[i]);
     337                        new_trans_input->AddTimeInput(xDynamicCast<TriaInput*>(input),transient_input->timesteps[i]);
     338                }
     339        }
     340        this->values=new_trans_input;
     341
     342        //      this->values->Echo();
     343        //this->values->Echo();
     344        //new_trans_input->Echo();
    254345
    255346}/*}}}*/

  • issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h

    r22519 r22739  
    1818                int    control_id;
    1919                int    enum_type;
     20                int      layout_enum;
    2021                Input* gradient;
    2122                Input* maxvalues;
    … …  
    2728                ControlInput();
    2829                ControlInput(int enum_type,int enum_input,IssmDouble* pvalues,IssmDouble* pmin,IssmDouble* pmax,int interp,int id);
     30                ControlInput(int enum_type,int enum_input,Input* pvalues,Input* pmin,Input* pmax,Input* pgrad,int interp,int id);
    2931                ~ControlInput();
    3032                /*}}}*/
    … …  
    7577                void Scale(IssmDouble scale_factor){_error_("not implemented yet");};
    7678                void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
     79                void SetGradient(Input* gradient_in,int timestep);
    7780                void SetGradient(Input* gradient_in);
    7881                void SetInput(Input* in_input);
     82                void SetInput(Input* in_input, int timeoffset);
    7983                void UpdateValue(IssmDouble scalar);
    8084                /*}}}*/

  • issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp

    r22519 r22739  
    8383
    8484        int i;
     85        printf("-------------- file: TransientInput.cpp line: %i\n",__LINE__);
    8586
    8687        _printf_("TransientInput:\n");
    … …  
    232233
    233234        /*First, recover current time from parameters: */
    234         this->parameters->FindParam(&time,TimeEnum);
     235        parameters->FindParam(&time,TimeEnum);
    235236
    236237        /*Retrieve interpolated values for this time step: */
    … …  
    268269/*}}}*/
    269270IssmDouble  TransientInput::GetTimeByOffset(int offset){/*{{{*/
    270 
    271271        if (offset < 0) offset=0;
    272         _assert_(offset<this->numtimesteps);
    273 
     272        _assert_(offset<(this->numtimesteps));
    274273        return this->timesteps[offset];
    275274}
    … …  
    469468        int        found;
    470469        int        offset;
    471         bool       interp;
     470        bool       interp=false;
    472471
    473472        /*First, recover interp bool: */

  • issm/trunk-jpl/src/c/cores/controladm1qn3_core.cpp

    r22633 r22739  
    6969        int num_responses,num_controls,numberofvertices,solution_type;
    7070        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    71        
     71        int* N = NULL;
     72        int N_add = 0;
     73
    7274        if (solution_type == TransientSolutionEnum){
    7375                femmodel = input_struct->femmodel->copy();
    74                 }
     76        }
    7577
    7678        IssmPDouble  *Jlist        = input_struct->Jlist;
    … …  
    8587        femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    8688        femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
     89        femmodel->parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
    8790        numberofvertices=femmodel->vertices->NumberOfVertices();
    8891
    … …  
    9295        GetPassiveVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
    9396        GetPassiveVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
    94         for(int i=0;i<numberofvertices;i++){
    95                 for(int c=0;c<num_controls;c++){
    96                         int index = num_controls*i+c;
     97
     98        N_add = 0;
     99        for (int c=0;c<num_controls;c++){
     100                for(int i=0;i<numberofvertices*N[c];i++){
     101                        int index = N_add*numberofvertices+i;
    97102                        X[index] = X[index]*reCast<double>(scaling_factors[c]);
    98103                        if(X[index]>XU[index]) X[index]=XU[index];
    … …  
    103108        /*Start Tracing*/
    104109        simul_starttrace(femmodel);
    105 
    106110        /*Set X as our new control input and as INDEPENDENT*/
    107111#ifdef _HAVE_AD_
    108         IssmDouble* aX=xNew<IssmDouble>(num_controls*numberofvertices,"t");
     112        IssmDouble* aX=xNew<IssmDouble>(intn,"t");
    109113#else
    110         IssmDouble* aX=xNew<IssmDouble>(num_controls*numberofvertices);
     114        IssmDouble* aX=xNew<IssmDouble>(intn);
    111115#endif
    112116        if(my_rank==0){
    … …  
    115119                }
    116120        }
     121
    117122        ISSM_MPI_Bcast(aX,intn,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    118123        SetControlInputsFromVectorx(femmodel,aX);
    … …  
    127132        if(solution_type==TransientSolutionEnum){
    128133                IssmDouble restart_time;
    129        
    130134                femmodel->parameters->FindParam(&restart_time,TimesteppingStartTimeEnum);
    131135                femmodel->parameters->SetParam(restart_time,TimeEnum);
    … …  
    139143        DataSet*    dependent_objects=NULL;
    140144        IssmDouble      J=0.;
    141        
    142145        femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
    143146        femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
    … …  
    262265
    263266                /*Add to totalgradient: */
    264                 if(my_rank==0) for(int i=0;i<num_independents;i++) totalgradient[i]+=weightVectorTimesJac[i];
     267                if(my_rank==0) for(int i=0;i<num_independents;i++) {
     268                        totalgradient[i]+=weightVectorTimesJac[i];
     269                }
    265270
    266271                /*free resources :*/
    … …  
    299304
    300305        /*Compute gradient*/
    301         for(long i=0;i<num_independents_old;i++)        G[i] = totalgradient[i];
     306        for(long i=0;i<num_independents_old;i++) G[i] = totalgradient[i];
    302307
    303308        /*Constrain Gradient*/
    304309        IssmDouble  Gnorm = 0.;
    305         for(int i=0;i<numberofvertices;i++){
    306                 for(int c=0;c<num_controls;c++){
    307                         int index = num_controls*i+c;
     310        N_add = 0;
     311        for (int c=0;c<num_controls;c++){
     312                for(int i=0;i<numberofvertices*N[c];i++){
     313                        int index = N_add*numberofvertices+i;
    308314                        if(X[index]>=XU[index]) G[index]=0.;
    309315                        if(X[index]<=XL[index]) G[index]=0.;
    … …  
    339345        double      *X  = NULL;
    340346        double      *G  = NULL;
     347        int*                    N       = NULL;
     348        int offset = 0;
     349        int N_add;
    341350
    342351        /*Recover some parameters*/
    … …  
    350359        femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
    351360        femmodel->parameters->SetParam(false,SaveResultsEnum);
     361        femmodel->parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
    352362        numberofvertices=femmodel->vertices->NumberOfVertices();
    353363
    … …  
    377387        Xpetsc->GetSize(&intn);
    378388        delete Xpetsc;
    379         _assert_(intn==numberofvertices*num_controls);
     389        //_assert_(intn==numberofvertices*num_controls);
    380390 
    381391        /*Get problem dimension and initialize gradient and initial guess*/
    … …  
    384394
    385395        /*Scale control for M1QN3*/
    386         for(int i=0;i<numberofvertices;i++){
    387                 for(int c=0;c<num_controls;c++){
    388                         int index = num_controls*i+c;
     396        N_add = 0;
     397        for (int c=0;c<num_controls;c++){
     398                for(int i=0;i<numberofvertices*N[c];i++){
     399                        int index = N_add*numberofvertices+i;
    389400                        X[index] = X[index]/reCast<IssmPDouble>(scaling_factors[c]);
    390401                }
    391         }
    392 
     402                N_add+=N[c];
     403        }
     404       
    393405        /*Allocate m1qn3 working arrays (see documentation)*/
    394406        long      m   = 100;
    … …  
    437449        GetPassiveVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
    438450
    439         for(int i=0;i<numberofvertices;i++){
    440                 for(int c=0;c<num_controls;c++){
    441                         int index = num_controls*i+c;
     451        N_add = 0;
     452        for (int c=0;c<num_controls;c++){
     453                for(int i=0;i<numberofvertices*N[c];i++){
     454                        int index = N_add*numberofvertices+i;
    442455                        X[index] = X[index]*reCast<double>(scaling_factors[c]);
    443456                        if(X[index]>XU[index]) X[index]=XU[index];
    444457                        if(X[index]<XL[index]) X[index]=XL[index];
    445458                }
     459                N_add +=N[c];
    446460        }
    447461               
    448462        /*Set X as our new control*/
    449463        IssmDouble* aX=xNew<IssmDouble>(intn);
    450         for(int i=0;i<intn;i++) aX[i] = reCast<IssmDouble>(X[i]);
     464        for(int i=0;i<intn;i++) {
     465                aX[i] = reCast<IssmDouble>(X[i]);
     466                }
    451467        SetControlInputsFromVectorx(femmodel,aX);
    452468        xDelete(aX);
    … …  
    454470        /*Set final gradient in inputs*/
    455471        IssmDouble* aG=xNew<IssmDouble>(intn);
    456         for(int i=0;i<intn;i++) aG[i] = reCast<IssmDouble>(G[i]);
     472        for(int i=0;i<intn;i++) {
     473                aG[i] = reCast<IssmDouble>(G[i]);
     474        }
    457475        ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG);
    458         xDelete(aG);
    459 
    460         /*Add last cost function to results*/
    461476
    462477        if (solution_type == TransientSolutionEnum){
    … …  
    465480                femmodel->OutputControlsx(&femmodel->results);
    466481                femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,mystruct.Jlist,(*mystruct.i),mystruct.N,1,0));
     482                femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,Gradient1Enum,G,numberofvertices,intn/numberofvertices,1,0));
     483                femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,InversionControlParametersEnum,X,numberofvertices,intn/numberofvertices,1,0));
    467484        }
    468485        else{
    … …  
    471488        femmodel->OutputControlsx(&femmodel->results);
    472489        }
     490
     491        xDelete(aG);
     492
     493        /*Add last cost function to results*/
     494
    473495        /*Finalize*/
    474496        if(VerboseControl()) _printf0_("   preparing final solution\n");

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

    r14999 r22739  
    99void ControlInputSetGradientx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,IssmDouble* gradient){
    1010
    11         /*Intermediaries*/
    12         int  num_controls;
    13         int *control_type = NULL;
     11        bool isautodiff;
     12        parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
     13        if(isautodiff){
     14                /*Intermediaries*/
     15                int  num_controls;
     16                int *control_type = NULL;
     17                int* M_all;
     18                int* N_all;
    1419
    15         /*Retrieve some parameters*/
    16         parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    17         parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
     20                /*Retrieve some parameters*/
     21                parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     22                parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
     23                parameters->FindParam(&M_all,NULL,ControlInputSizeMEnum);
     24                parameters->FindParam(&N_all,NULL,ControlInputSizeNEnum);
    1825
    19         for(int i=0;i<num_controls;i++){
    20                 for(int j=0;j<elements->Size();j++){
    21                         Element* element=(Element*)elements->GetObjectByOffset(j);
    22                         element->ControlInputSetGradient(gradient,control_type[i],i);
     26                int offset = 0;
     27                for(int i=0;i<num_controls;i++){
     28                        for(int j=0;j<elements->Size();j++){
     29                                Element* element=(Element*)elements->GetObjectByOffset(j);
     30                                element->ControlInputSetGradient(gradient,control_type[i],i,offset,N_all[i],M_all[i]);
     31                        }
     32                        offset+=M_all[i]*N_all[i];
    2333                }
     34
     35                /*Clean up and return*/
     36                xDelete<int>(control_type);
    2437        }
     38        else{
     39                int  num_controls;
     40                int *control_type = NULL;
    2541
    26         /*Clean up and return*/
    27         xDelete<int>(control_type);
     42                /*Retrieve some parameters*/
     43                parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     44                parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
     45
     46                int offset = 0;
     47                for(int i=0;i<num_controls;i++){
     48                        for(int j=0;j<elements->Size();j++){
     49                                Element* element=(Element*)elements->GetObjectByOffset(j);
     50                                element->ControlInputSetGradient(gradient,control_type[i],i);
     51                        }
     52                }
     53
     54                /*Clean up and return*/
     55                xDelete<int>(control_type);
     56        }
    2857
    2958}

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

    r22424 r22739  
    99void GetVectorFromControlInputsx(Vector<IssmDouble>** pvector, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,const char* data,bool onsid){/*{{{*/
    1010
    11         int  num_controls;
    12         int *control_type = NULL;
    13         Vector<IssmDouble>*  vector=NULL;
     11        bool isautodiff;
     12        parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
     13        if(isautodiff){
     14                int*  N = NULL;
     15                int*  M = NULL;
     16                int  num_controls;
     17                int* control_type = NULL;
     18                Vector<IssmDouble>*  vector=NULL;
    1419
    15         /*Retrieve some parameters*/
    16         parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    17         parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
     20                /*Retrieve some parameters*/
     21                parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     22                parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
     23                parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
     24                parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
    1825
    19         /*Allocate and populate gradient*/
    20         vector=new Vector<IssmDouble>(num_controls*vertices->NumberOfVertices());
     26                /*1. Get vector size*/
     27                int size = 0;
     28                for(int i=0;i<num_controls;i++) size+=M[i]*N[i];
    2129
    22         for(int i=0;i<num_controls;i++){
    23                 for(int j=0;j<elements->Size();j++){
    24                         Element* element=(Element*)elements->GetObjectByOffset(j);
    25                         element->GetVectorFromControlInputs(vector,control_type[i],i,data,onsid);
     30
     31                /*2. Allocate vector*/
     32                vector=new Vector<IssmDouble>(size);
     33
     34                /*3. Populate vector*/
     35                int offset = 0;
     36                for(int i=0;i<num_controls;i++){
     37                        for(int j=0;j<elements->Size();j++){
     38                                Element* element=(Element*)elements->GetObjectByOffset(j);
     39                                element->GetVectorFromControlInputs(vector,control_type[i],i,data,offset,onsid);
     40                        }
     41                        offset += M[i]*N[i];
    2642                }
     43
     44                vector->Assemble();
     45
     46                /*Assign output pointers:*/
     47                xDelete<int>(control_type);
     48                xDelete<int>(M);
     49                xDelete<int>(N);
     50
     51                *pvector=vector;
     52        }
     53        else{
     54                int  num_controls;
     55                int* control_type = NULL;
     56                Vector<IssmDouble>*  vector=NULL;
     57
     58                /*Retrieve some parameters*/
     59                parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     60                parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
     61
     62
     63                /*2. Allocate vector*/
     64                vector=new Vector<IssmDouble>(num_controls*vertices->NumberOfVertices());
     65
     66                /*3. Populate vector*/
     67                int offset = 0;
     68                for(int i=0;i<num_controls;i++){
     69                        for(int j=0;j<elements->Size();j++){
     70                                Element* element=(Element*)elements->GetObjectByOffset(j);
     71                                element->GetVectorFromControlInputs(vector,control_type[i],i,data,onsid);
     72                        }
     73                }
     74                vector->Assemble();
     75
     76                /*Assign output pointers:*/
     77                xDelete<int>(control_type);
     78                *pvector=vector;
    2779        }
    2880
    29         vector->Assemble();
    30 
    31         /*Assign output pointers:*/
    32         xDelete<int>(control_type);
    33         *pvector=vector;
    3481}/*}}}*/
    3582void GetVectorFromControlInputsx( IssmDouble** pvector, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters, const char* data,bool onsid){/*{{{*/
    36 
    37         /*output: */
    38         IssmDouble* vector=NULL;
    3983
    4084        /*intermediary: */
    … …  
    4286
    4387        GetVectorFromControlInputsx( &vec_vector, elements,nodes, vertices, loads, materials, parameters,data,onsid);
    44         vector=vec_vector->ToMPISerial();
     88        IssmDouble* vector=vec_vector->ToMPISerial();
    4589
    4690        /*Free ressources:*/

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp

    r22536 r22739  
    1010
    1111
    12 #if !defined(_HAVE_ADOLC_)
    13 void    UpdateElementsAndMaterialsControl(Elements* elements,Materials* materials, IoModel* iomodel){
    14        
     12//#if !defined(_HAVE_ADOLC_)
     13void    UpdateElementsAndMaterialsControl(Elements* elements,Parameters* parameters,Materials* materials, IoModel* iomodel){
     14        bool isautodiff;
     15        iomodel->FindConstant(&isautodiff,"md.autodiff.isautodiff");
     16        if(!isautodiff){
    1517        /*Intermediary*/
    1618        bool       control_analysis;
    … …  
    6567                }
    6668        }
     69        parameters->AddObject(new IntParam(ControlInputSizeMEnum,iomodel->numberofvertices));
    6770
    6871        for(int i=0;i<num_controls;i++){
    … …  
    128131        xDelete<char*>(controls);
    129132}
    130 #else
    131 void UpdateElementsAndMaterialsControl(Elements* elements,Materials* materials, IoModel* iomodel){
    132 
     133//}
     134//#else
     135//void UpdateElementsAndMaterialsControl(Elements* elements,Parameters* parameters,Materials* materials, IoModel* iomodel){
     136else{
    133137
    134138        /*Intermediaries*/
    135         int                             num_independent_objects,M,N;
     139        int                             num_independent_objects,M,N,M_par,N_par;
    136140        char**                  names                   = NULL;
    137141        int*                            types                                                   = NULL;
     142        int*                            control_sizes                           = NULL;
     143        int*                            M_all                                                   = NULL;
    138144        IssmDouble*             independent                                     = NULL;
    139         IssmDouble**    independents_min                        = NULL;
    140         IssmDouble**    independents_max                        = NULL;
     145        IssmDouble*             independents_fullmin    = NULL;
     146        IssmDouble*             independents_fullmax            = NULL;
    141147        bool                            control_analysis                        =false;
    142148
    143149        iomodel->FindConstant(&control_analysis,"md.inversion.iscontrol");
    144        
     150
    145151        /*Now, return if no control*/
    146152        if(!control_analysis) return;
    … …  
    154160
    155161               
    156         /*TODO: fetch min and max*/
    157         independents_min = xNew<IssmDouble*>(num_independent_objects);
    158         independents_max = xNew<IssmDouble*>(num_independent_objects);
    159         for(int i=0;i<num_independent_objects;i++) independents_min[i]=NULL;
    160         for(int i=0;i<num_independent_objects;i++) independents_max[i]=NULL;
     162        M_all = xNew<int>(num_independent_objects);
    161163
    162164        /*create independent objects, and at the same time, fetch the corresponding independent variables,
    163165         *and declare them as such in ADOLC: */
     166        iomodel->FetchData(&independents_fullmin,&M_par,&N_par,"md.autodiff.independent_min_parameters");
     167        iomodel->FetchData(&independents_fullmax,&M_par,&N_par,"md.autodiff.independent_max_parameters");
     168        iomodel->FetchData(&control_sizes,NULL,NULL,"md.autodiff.independent_control_sizes");
     169       
     170        int* start_point = NULL;
     171        start_point = xNew<int>(num_independent_objects);
     172        int counter = 0;
     173        for(int i=0;i<num_independent_objects;i++){
     174                start_point[i]=counter;
     175                counter+=control_sizes[i];
     176        }
     177
    164178        for(int i=0;i<num_independent_objects;i++){
    165179
    … …  
    169183                        char* iofieldname  = NULL;
    170184                        int   input_enum;
     185                        IssmDouble*     independents_min                        = NULL;
     186                        IssmDouble*        independents_max                     = NULL;
     187       
    171188                        FieldAndEnumFromCode(&input_enum,&iofieldname,names[i]);
    172189
    173190                        /*Fetch required data*/
    174191                        iomodel->FetchData(&independent,&M,&N,iofieldname);
    175                         iomodel->FetchData(&independents_min[i],&M,&N,"md.autodiff.independent_min_parameters");
    176                         iomodel->FetchData(&independents_max[i],&M,&N,"md.autodiff.independent_max_parameters");
    177                        
    178                         _assert_(independent);
     192                        _assert_(independent);
     193                        _assert_(N==control_sizes[i]);
     194                        _printf_("N control size: "<<control_sizes[i]<<"\n");
     195               
     196                        independents_min = NULL; independents_min = xNew<IssmDouble>(M*N);
     197                        independents_max = NULL; independents_max = xNew<IssmDouble>(M*N);
     198                        for(int m=0;m<M;m++){
     199                                for(int n=0;n<N;n++){
     200                                        independents_min[N*m+n]=independents_fullmin[N_par*m+start_point[i]+n];
     201                                        independents_max[N*m+n]=independents_fullmax[N_par*m+start_point[i]+n];
     202                                }
     203                        }
     204                        if(N!=1) M_all[i]=M-1;
     205                        else M_all[i]=M;
    179206
    180207                        for(int j=0;j<elements->Size();j++){
    181208                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
    182                                 element->ControlInputCreate(independent,independents_min[i],independents_max[i],iomodel,M,N,input_enum,i+1);
     209                                element->ControlInputCreate(independent,independents_min,independents_max,iomodel,M,N,input_enum,i+1);
    183210                        }
    184211                        xDelete<IssmDouble>(independent);
     212                        xDelete<IssmDouble>(independents_min);
     213                        xDelete<IssmDouble>(independents_max);
     214       
    185215                }
    186216                else{
    … …  
    188218                }
    189219        }
     220        parameters->AddObject(new IntVecParam(ControlInputSizeNEnum,control_sizes,num_independent_objects));
     221        parameters->AddObject(new IntVecParam(ControlInputSizeMEnum,M_all,num_independent_objects));
    190222
    191223        /*cleanup*/
    192224        for(int i=0;i<num_independent_objects;i++){
    193225                xDelete<char>(names[i]);
    194                 xDelete<IssmDouble>(independents_min[i]);
    195                 xDelete<IssmDouble>(independents_max[i]);
    196         }
    197         xDelete<IssmDouble*>(independents_min);
    198         xDelete<IssmDouble*>(independents_max);
     226        }
    199227        xDelete<char*>(names);
    200228        xDelete<int>(types);
    201 
     229        xDelete<int>(M_all);
     230        xDelete<IssmDouble>(independents_fullmin);
     231        xDelete<IssmDouble>(independents_fullmax);
     232        xDelete<int>(start_point);
    202233        /*Step2: create cost functions (dependents)*/
    203234
    204235        return;
    205236}
    206 #endif
     237}
     238//#endif

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

    r22004 r22739  
    6969        /*Solution specific updates*/
    7070        if(VerboseMProcessor()) _printf0_("   updating elements and materials for control parameters" << "\n");
    71         UpdateElementsAndMaterialsControl(elements,materials,iomodel);
     71        UpdateElementsAndMaterialsControl(elements,parameters,materials,iomodel);
    7272        #ifdef _HAVE_DAKOTA_
    7373        if(VerboseMProcessor()) _printf0_("   updating elements and materials for uncertainty quantification" << "\n");

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

    r22004 r22739  
    1818void CreateParametersDakota(Parameters* parameters,IoModel* iomodel,char* rootpath);
    1919void CreateOutputDefinitions(Elements* elements, Parameters* parameters,IoModel* iomodel);
    20 void UpdateElementsAndMaterialsControl(Elements* elements,Materials* materials, IoModel* iomodel);
     20void UpdateElementsAndMaterialsControl(Elements* elements,Parameters* parameters,Materials* materials, IoModel* iomodel);
    2121void UpdateElementsAndMaterialsDakota(Elements* elements,Materials* materials, IoModel* iomodel);
    2222void UpdateElementsTransient(Elements* elements,Parameters* parameters,IoModel* iomodel,int analysis_type);

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

    r18128 r22739  
    99void SetControlInputsFromVectorx(FemModel* femmodel,IssmDouble* vector){
    1010
    11         int  num_controls;
    12         int *control_type = NULL;
     11        bool isautodiff;
     12        femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
     13        if(isautodiff){
     14                int  num_controls;
     15                int* control_type = NULL;
     16                int* M = NULL;
     17                int* N = NULL;
    1318
    14         /*Retrieve some parameters*/
    15         femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    16         femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
     19                /*Retrieve some parameters*/
     20                femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     21                femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
     22                femmodel->parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
     23                femmodel->parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
    1724
    18         for(int i=0;i<num_controls;i++){
    19                 for(int j=0;j<femmodel->elements->Size();j++){
    20                         Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
    21                         element->SetControlInputsFromVector(vector,control_type[i],i);
     25                int offset = 0;
     26                for(int i=0;i<num_controls;i++){
     27                        for(int j=0;j<femmodel->elements->Size();j++){
     28                                Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
     29                                element->SetControlInputsFromVector(vector,control_type[i],i,offset,N[i],M[i]);
     30                        }
     31                offset += M[i]*N[i];
    2232                }
     33
     34
     35                xDelete<int>(control_type);
    2336        }
     37        else{
    2438
    25         xDelete<int>(control_type);
     39                int  num_controls;
     40                int* control_type = NULL;
     41                femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     42                femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
     43                int offset = 0;
     44                for(int i=0;i<num_controls;i++){
     45                        for(int j=0;j<femmodel->elements->Size();j++){
     46                                Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
     47                                element->SetControlInputsFromVector(vector,control_type[i],i);
     48                        }
     49                }
     50                xDelete<int>(control_type);
     51        }
    2652}
    2753
Note: See TracChangeset for help on using the changeset viewer.