Changeset 22515


Ignore:
Timestamp:
03/10/18 19:14:51 (7 years ago)
Author:
erobo
Message:

CHG: changes to allow AD/m1qn3 to work

Location:
issm/trunk-jpl/src
Files:
19 edited

Legend:

Unmodified
Added
Removed

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

    r22507 r22515  
    119119                 IssmDouble J=0.;
    120120                 IssmDouble J_sum=0.;
    121        
     121        printf("-------------- file: Cfsurfacesquare.cpp line: %i\n",__LINE__);
     122        _printf_("current time: "<<time<<"   datatime: "<<datatime<<"\n");
    122123         if(datatime<=time && !timepassedflag){
     124printf("-------------- file: Cfsurfacesquare.cpp line: %i\n",__LINE__);
    123125                 for(i=0;i<femmodel->elements->Size();i++){
    124126                         Element* element=(Element*)femmodel->elements->GetObjectByOffset(i);

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

    r22507 r22515  
    16021602    }
    16031603    else _error_("Cannot add input for vector type " << vector_type << " (not supported)");
     1604}
     1605/*}}}*/
     1606void       Element::ControlInputCreate(IssmDouble* vector,IssmDouble* min_vector,IssmDouble* max_vector,IoModel* iomodel,int M,int N,int input_enum,int id){/*{{{*/
     1607
     1608        /*Intermediaries*/
     1609        int         numvertices = this->GetNumberOfVertices();
     1610        int        *vertexids   = xNew<int>(numvertices);
     1611        IssmDouble *values      = xNew<IssmDouble>(numvertices);
     1612        IssmDouble *values_min  = xNew<IssmDouble>(numvertices);
     1613        IssmDouble *values_max  = xNew<IssmDouble>(numvertices);
     1614
     1615        /*Some sanity checks*/
     1616        _assert_(vector);
     1617        _assert_(min_vector);
     1618        _assert_(max_vector);
     1619
     1620        /*For now we only support nodal vectors*/
     1621        if(M!=iomodel->numberofvertices) _error_("not supported");
     1622        if(N!=1) _error_("not supported");
     1623
     1624        /*Recover vertices ids needed to initialize inputs*/
     1625        _assert_(iomodel->elements);
     1626        for(int i=0;i<numvertices;i++){
     1627                vertexids[i]=reCast<int>(iomodel->elements[numvertices*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
     1628        }
     1629
     1630        /*Are we in transient or static? */
     1631        if(M==iomodel->numberofvertices){
     1632                for(int i=0;i<numvertices;i++){
     1633                        values[i]=vector[vertexids[i]-1];
     1634                        values_min[i] = min_vector[vertexids[i]-1];
     1635                        values_max[i] = max_vector[vertexids[i]-1];
     1636                }
     1637                this->AddControlInput(input_enum,values,min_vector,max_vector,P1Enum,id);
     1638        }
     1639
     1640        /*clean up*/
     1641        xDelete<IssmDouble>(values);
     1642        xDelete<int>(vertexids);
    16041643}
    16051644/*}}}*/

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

    r22507 r22515  
    122122                void               InputChangeName(int enum_type,int enum_type_old);
    123123                void               InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code);
     124                void               ControlInputCreate(IssmDouble* doublearray,IssmDouble* independents_min,IssmDouble* independents_max,IoModel* iomodel,int M,int N,int input_enum,int id);
    124125                void                                     DatasetInputAdd(int enum_type,IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code,int input_enum);
    125126                void               InputDuplicate(int original_enum,int new_enum);
    … …  
    188189                virtual void       AddBasalInput(int input_enum, IssmDouble* values, int interpolation_enum)=0;
    189190                virtual void       AddInput(int input_enum, IssmDouble* values, int interpolation_enum)=0;
     191                virtual void       AddControlInput(int input_enum, IssmDouble* values,IssmDouble* values_min,IssmDouble* values_max, int interpolation_enum,int id){_error_("not supported yet");};
     192                virtual void       DatasetInputCreate(IssmDouble* array,int M,int N,int* individual_enums,int num_inputs,IoModel* iomodel,int input_enum){_error_("not supported");};
    190193                virtual void       AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part)=0;
    191194                virtual void             BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement){_error_("not implemented yet");};

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

    r22379 r22515  
    14001400                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
    14011401                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
    1402                                                 this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1402                                                this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1));
    14031403                                        }
    14041404                                        break;
    … …  
    14081408                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
    14091409                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
    1410                                                 this->inputs->AddInput(new ControlInput(VxEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1410                                                this->inputs->AddInput(new ControlInput(VxEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1));
    14111411                                        }
    14121412                                        break;
    … …  
    14161416                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
    14171417                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
    1418                                                 this->inputs->AddInput(new ControlInput(VyEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1418                                                this->inputs->AddInput(new ControlInput(VyEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1));
    14191419                                        }
    14201420                                        break;
    … …  
    14241424                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i];
    14251425                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i];
    1426                                                 this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1426                                                this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1));
    14271427                                        }
    14281428                                        break;
    … …  
    14321432                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i];
    14331433                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i];
    1434                                                 this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1434                                                this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1));
    14351435                                        }
    14361436                                        break;
    … …  
    14401440                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i];
    14411441                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i];
    1442                                                 this->inputs->AddInput(new ControlInput(DamageDEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1442                                                this->inputs->AddInput(new ControlInput(DamageDEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1));
    14431443                                        }
    14441444                                        break;

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

    r21876 r22515  
    394394                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts;
    395395                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts;
    396                                                 this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     396                                                this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1));
    397397                                        }
    398398                                        break;
    … …  
    402402                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts;
    403403                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts;
    404                                                 this->inputs->AddInput(new ControlInput(VxEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     404                                                this->inputs->AddInput(new ControlInput(VxEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1));
    405405                                        }
    406406                                        break;
    … …  
    410410                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts;
    411411                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i]/yts;
    412                                                 this->inputs->AddInput(new ControlInput(VyEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     412                                                this->inputs->AddInput(new ControlInput(VyEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1));
    413413                                        }
    414414                                        break;
    … …  
    418418                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i];
    419419                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i];
    420                                                 this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     420                                                this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1));
    421421                                        }
    422422                                        break;
    … …  
    426426                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i];
    427427                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i];
    428                                                 this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     428                                                this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1));
    429429                                        }
    430430                                        break;
    … …  
    434434                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i];
    435435                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(tetra_vertex_ids[j]-1)*num_control_type+i];
    436                                                 this->inputs->AddInput(new ControlInput(DamageDEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     436                                                this->inputs->AddInput(new ControlInput(DamageDEnum,TetraInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1));
    437437                                        }
    438438                                        break;

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

    r22488 r22515  
    168168        _assert_(this->inputs);
    169169        this->inputs->AddInput(new TriaInput(input_enum,values,interpolation_enum));
     170}
     171/*}}}*/
     172void       Tria::AddControlInput(int input_enum,IssmDouble* values,IssmDouble* values_min,IssmDouble* values_max, int interpolation_enum,int id){/*{{{*/
     173
     174        /*Call inputs method*/
     175        _assert_(this->inputs);
     176        this->inputs->AddInput(new ControlInput(input_enum,TriaInputEnum,values,values_min,values_max,interpolation_enum,id));
     177}
     178/*}}}*/
     179void       Tria::DatasetInputCreate(IssmDouble* array,int M,int N,int* individual_enums,int num_inputs,IoModel* iomodel,int input_enum){/*{{{*/
     180
     181        IssmDouble nodeinputs[NUMVERTICES];
     182        if(num_inputs<1) _error_("Cannot create a DatasetInput of size <1");
     183        if(M!=iomodel->numberofvertices) _error_("not supported yet");
     184        if(N!=num_inputs) _error_("sizes are not consistent");
     185
     186        int        tria_vertex_ids[3];
     187       
     188        for(int k=0;k<3;k++){
     189                tria_vertex_ids[k]=reCast<int>(iomodel->elements[3*this->Sid()+k]); //ids for vertices are in the elements array from Matlab
     190        }
     191        /*Create inputs and add to DataSetInput*/
     192        DatasetInput* datasetinput=new DatasetInput(input_enum);
     193        for(int i=0;i<num_inputs;i++){
     194                for(int j=0;j<NUMVERTICES;j++)nodeinputs[j]=array[(tria_vertex_ids[j]-1)*N+i];
     195                datasetinput->AddInput(new TriaInput(input_enum,nodeinputs,P1Enum),individual_enums[i]);
     196        }
     197
     198        /*Add datasetinput to element inputs*/
     199        this->inputs->AddInput(datasetinput);
    170200}
    171201/*}}}*/
    … …  
    19111941        IssmDouble cmmininputs[3];
    19121942        IssmDouble cmmaxinputs[3];
    1913         bool       control_analysis   = false;
     1943        bool       control_analysis,ad_analysis   = false;
    19141944        int        num_control_type,num_responses;
    19151945        char**     controls = NULL;
    … …  
    19191949        iomodel->FindConstant(&yts,"md.constants.yts");
    19201950        iomodel->FindConstant(&control_analysis,"md.inversion.iscontrol");
    1921         if(control_analysis) iomodel->FindConstant(&num_control_type,"md.inversion.num_control_parameters");
    1922         if(control_analysis) iomodel->FindConstant(&num_responses,"md.inversion.num_cost_functions");
     1951        iomodel->FindConstant(&ad_analysis, "md.autodiff.isautodiff");
     1952        if(control_analysis && !ad_analysis) iomodel->FindConstant(&num_control_type,"md.inversion.num_control_parameters");
     1953        if(control_analysis && !ad_analysis) iomodel->FindConstant(&num_responses,"md.inversion.num_cost_functions");
     1954        if(control_analysis && ad_analysis) iomodel->FindConstant(&num_control_type,"md.autodiff.num_independent_objects");
     1955        if(control_analysis && ad_analysis) iomodel->FindConstant(&num_responses,"md.autodiff.num_dependent_objects");
     1956
     1957
    19231958
    19241959        /*Recover vertices ids needed to initialize inputs*/
    … …  
    19341969        /*Control Inputs*/
    19351970        if (control_analysis){
    1936                 iomodel->FindConstant(&controls,NULL,"md.inversion.control_parameters");
     1971                if(!ad_analysis)iomodel->FindConstant(&controls,NULL,"md.inversion.control_parameters");
     1972                if(ad_analysis)iomodel->FindConstant(&controls,NULL,"md.autodiff.independent_object_names");
     1973
    19371974                for(i=0;i<num_control_type;i++){
    19381975                        _assert_(controls[i]);
    … …  
    19441981                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
    19451982                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
    1946                                                 this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1983                                                this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1));
    19471984                                        }
    19481985                                        break;
    … …  
    19521989                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
    19531990                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
    1954                                                 this->inputs->AddInput(new ControlInput(VxEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1991                                                this->inputs->AddInput(new ControlInput(VxEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1));
    19551992                                        }
    19561993                                        break;
    … …  
    19601997                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
    19611998                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
    1962                                                 this->inputs->AddInput(new ControlInput(VyEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1999                                                this->inputs->AddInput(new ControlInput(VyEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1));
    19632000                                        }
    19642001                                        break;
    … …  
    19682005                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i];
    19692006                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i];
    1970                                                 this->inputs->AddInput(new ControlInput(ThicknessEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     2007                                                this->inputs->AddInput(new ControlInput(ThicknessEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1));
    19712008                                        }
    19722009                                        break;
    … …  
    19762013                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i];
    19772014                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i];
    1978                                                 this->inputs->AddInput(new ControlInput(BalancethicknessSpcthicknessEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     2015                                                this->inputs->AddInput(new ControlInput(BalancethicknessSpcthicknessEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1));
    19792016                                        }
    19802017                                        break;
    … …  
    19842021                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i];
    19852022                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i];
    1986                                                 this->inputs->AddInput(new ControlInput(BalancethicknessOmegaEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     2023                                                this->inputs->AddInput(new ControlInput(BalancethicknessOmegaEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1));
    19872024                                        }
    19882025                                        break;
    … …  
    19922029                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i];
    19932030                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i];
    1994                                                 this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     2031                                                this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1));
    19952032                                        }
    19962033                                        break;
    … …  
    20002037                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i];
    20012038                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i];
    2002                                                 this->inputs->AddInput(new ControlInput(MaterialsRheologyBbarEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     2039                                                this->inputs->AddInput(new ControlInput(MaterialsRheologyBbarEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1));
    20032040                                        }
    20042041                                        break;
    … …  
    20082045                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i];
    20092046                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i];
    2010                                                 this->inputs->AddInput(new ControlInput(DamageDbarEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     2047                                                this->inputs->AddInput(new ControlInput(DamageDbarEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,P1Enum,i+1));
    20112048                                        }
    20122049                                        break;

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

    r22473 r22515  
    162162                void           AddBasalInput(int input_enum, IssmDouble* values, int interpolation_enum);
    163163                void           AddInput(int input_enum, IssmDouble* values, int interpolation_enum);
     164                void           AddControlInput(int input_enum, IssmDouble* values,IssmDouble* values_min,IssmDouble* values_max, int interpolation_enum,int id);
     165                void           DatasetInputCreate(IssmDouble* array,int M,int N,int* individual_enums,int num_inputs,IoModel* iomodel,int input_enum);
    164166                IssmDouble     GetArea(void);
    165167                IssmDouble         GetArea3D(void);

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

    r21974 r22515  
    2222}
    2323/*}}}*/
    24 ControlInput::ControlInput(int in_enum_type,int enum_input,IssmDouble* pvalues,IssmDouble* pmin,IssmDouble* pmax,int id){/*{{{*/
     24ControlInput::ControlInput(int in_enum_type,int input_layout_enum,IssmDouble* pvalues,IssmDouble* pmin,IssmDouble* pmax,int interp,int id){/*{{{*/
    2525
    2626        control_id=id;
    2727        enum_type=in_enum_type;
    2828
    29         switch(enum_input){
     29        _assert_(interp==P1Enum);
     30
     31        switch(input_layout_enum){
    3032                case TriaInputEnum:
    31                         values     =new TriaInput(enum_type,pvalues,P1Enum);
    32                         savedvalues=new TriaInput(enum_type,pvalues,P1Enum);
    33                         minvalues  =new TriaInput(enum_type,pmin,P1Enum);
    34                         maxvalues  =new TriaInput(enum_type,pmax,P1Enum);
     33                        values     =new TriaInput(enum_type,pvalues,interp);
     34                        savedvalues=new TriaInput(enum_type,pvalues,interp);
     35                        minvalues  =new TriaInput(enum_type,pmin,interp);
     36                        maxvalues  =new TriaInput(enum_type,pmax,interp);
    3537                        break;
    3638                case PentaInputEnum:
    37                         values     =new PentaInput(enum_type,pvalues,P1Enum);
    38                         savedvalues=new PentaInput(enum_type,pvalues,P1Enum);
    39                         minvalues  =new PentaInput(enum_type,pmin,P1Enum);
    40                         maxvalues  =new PentaInput(enum_type,pmax,P1Enum);
     39                        values     =new PentaInput(enum_type,pvalues,interp);
     40                        savedvalues=new PentaInput(enum_type,pvalues,interp);
     41                        minvalues  =new PentaInput(enum_type,pmin,interp);
     42                        maxvalues  =new PentaInput(enum_type,pmax,interp);
    4143                        break;
    4244                default:
    43                         _error_("Input of Enum " << EnumToStringx(enum_input) << " not supported yet by ControlInput");
     45                        _error_("Input of Enum " << EnumToStringx(input_layout_enum) << " not supported yet by ControlInput");
    4446        }
    4547        gradient   =NULL;

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

    r21974 r22515  
    2626                /*ControlInput constructors, destructors: {{{*/
    2727                ControlInput();
    28                 ControlInput(int enum_type,int enum_input,IssmDouble* pvalues,IssmDouble* pmin,IssmDouble* pmax,int id);
     28                ControlInput(int enum_type,int enum_input,IssmDouble* pvalues,IssmDouble* pmin,IssmDouble* pmax,int interp,int id);
    2929                ~ControlInput();
    3030                /*}}}*/

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

    r22445 r22515  
    355355        long niter = long(maxsteps); /*Maximum number of iterations*/
    356356        long nsim  = long(maxiter);/*Maximum number of function calls*/
    357 
     357       
    358358        /*Get initial guess*/
    359359        Vector<double> *Xpetsc = NULL;
     360
     361        /*THIS IS WHERE IT FAILS*/
    360362        GetPassiveVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
    361363        X = Xpetsc->ToMPISerial();
    … …  
    363365        delete Xpetsc;
    364366        _assert_(intn==numberofvertices*num_controls);
    365 
     367 
    366368        /*Get problem dimension and initialize gradient and initial guess*/
    367369        long n = long(intn);
    … …  
    395397        /*Initialize Gradient and cost function of M1QN3*/
    396398        indic = 4; /*gradient required*/
     399printf("-------------- file: controladm1qn3_core.cpp line: %i\n",__LINE__);
    397400        simul_ad(&indic,&n,X,&f,G,izs,rzs,(void*)&mystruct);
    398401        /*Estimation of the expected decrease in f during the first iteration*/

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

    r22443 r22515  
    2828        if(control_analysis){
    2929
    30                 /*How many controls and how many responses?*/
    31                 parameters->AddObject(iomodel->CopyConstantObject("md.inversion.num_control_parameters",InversionNumControlParametersEnum));
    32                 parameters->AddObject(iomodel->CopyConstantObject("md.inversion.num_cost_functions",InversionNumCostFunctionsEnum));
    33 
    3430                /*What solution type?*/
    3531                if(solution_type==SteadystateSolutionEnum){
    … …  
    4036                }
    4137
    42                 /*recover controls and convert to Enums*/
    43                 iomodel->FindConstant(&controls,&num_controls,"md.inversion.control_parameters");
    44                 if(num_controls<1) _error_("no controls found");
    45                 int* control_enums=xNew<int>(num_controls);
    46                 for(int i=0;i<num_controls;i++){
    47                         control_enums[i]=StringToEnumx(controls[i]);
    48                         xDelete<char>(controls[i]);
     38                switch(inversiontype){
     39                        _printf_("Inversiontype: "<<inversiontype<<"\n");
     40                          {
     41                        case 0:/*Brent Search*/
     42                        case 1:/*TAO*/
     43                        case 2:/*M1QN3*/
     44                        case 3:/*Validation*/
     45                                /*How many controls and how many responses?*/
     46                                parameters->AddObject(iomodel->CopyConstantObject("md.inversion.num_control_parameters",InversionNumControlParametersEnum));
     47                                parameters->AddObject(iomodel->CopyConstantObject("md.inversion.num_cost_functions",InversionNumCostFunctionsEnum));
     48
     49                                /*recover controls and convert to Enums*/
     50                                iomodel->FindConstant(&controls,&num_controls,"md.inversion.control_parameters");
     51                                if(num_controls<1) _error_("no controls found");
     52                                int* control_enums=xNew<int>(num_controls);
     53                                for(int i=0;i<num_controls;i++){
     54                                        control_enums[i]=StringToEnumx(controls[i]);
     55                                        xDelete<char>(controls[i]);
     56                                }
     57                                xDelete<char*>(controls);
     58                                parameters->AddObject(new IntVecParam(InversionControlParametersEnum,control_enums,num_controls));
     59
     60                                iomodel->FindConstant(&cm_responses,&num_costfunc,"md.inversion.cost_functions");
     61                                if(num_costfunc<1) _error_ ("no cost functions found");
     62                                int* costfunc_enums=xNew<int>(num_costfunc);
     63                                for(int i=0;i<num_costfunc;i++){
     64                                        costfunc_enums[i]=StringToEnumx(cm_responses[i]);
     65                                        xDelete<char>(cm_responses[i]);
     66                                }
     67                                xDelete<char*>(cm_responses);
     68                                parameters->AddObject(new IntVecParam(InversionCostFunctionsEnum,costfunc_enums,num_costfunc));
     69
     70                                xDelete<int>(control_enums);
     71                                xDelete<int>(costfunc_enums);
     72
     73                                break;
    4974                }
    50                 xDelete<char*>(controls);
    51                 parameters->AddObject(new IntVecParam(InversionControlParametersEnum,control_enums,num_controls));
     75                        case 4:/*AD M1QN3*/
     76                        {
     77                        /*Intermediaries*/
     78                        int            num_independent_objects,M;
     79                        char**         names                   = NULL;
     80                               
     81                                /*this is done somewhere else*/
     82                                parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.num_independent_objects",InversionNumControlParametersEnum));
     83                           parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.num_dependent_objects",InversionNumCostFunctionsEnum));
     84                               
     85                                /*Step1: create controls (independents)*/
     86                                iomodel->FetchData(&num_independent_objects,"md.autodiff.num_independent_objects");
     87                                _assert_(num_independent_objects>0);
     88                                iomodel->FetchData(&names,&M,"md.autodiff.independent_object_names");
     89                                _assert_(M==num_independent_objects);
     90                                int* ind_enums=xNew<int>(num_independent_objects);
     91                                for(int i=0;i<num_independent_objects;i++){
     92                                        ind_enums[i]=StringToEnumx(names[i]);
     93                                        xDelete<char>(names[i]);
     94                                }
    5295
    53                 iomodel->FindConstant(&cm_responses,&num_costfunc,"md.inversion.cost_functions");
    54                 if(num_costfunc<1) _error_ ("no cost functions found");
    55                 int* costfunc_enums=xNew<int>(num_costfunc);
    56                 for(int i=0;i<num_costfunc;i++){
    57                         costfunc_enums[i]=StringToEnumx(cm_responses[i]);
    58                         xDelete<char>(cm_responses[i]);
     96                                parameters->AddObject(new IntVecParam(InversionControlParametersEnum,ind_enums,num_independent_objects));
     97                                //iomodel->FetchData(&num_costfunc,"md.numberedcostfunction.num_cost_functions");
     98                                      _assert_(num_costfunc>0);
     99                                iomodel->FindConstant(&cm_responses,&num_costfunc,"md.autodiff.dependent_object_names");
     100                                if(num_costfunc<1) _error_ ("no cost functions found");
     101                                int* costfunc_enums=xNew<int>(num_costfunc);
     102                                for(int i=0;i<num_costfunc;i++){
     103                                        costfunc_enums[i]=StringToEnumx(cm_responses[i]);
     104                                        xDelete<char>(cm_responses[i]);
     105                                }
     106                                xDelete<char*>(cm_responses);
     107                                parameters->AddObject(new IntVecParam(InversionCostFunctionsEnum,costfunc_enums,num_costfunc));
     108                               
     109                                iomodel->FetchData(&control_scaling_factors,NULL,NULL,"md.autodiff.independent_scaling_factors");
     110                                parameters->AddObject(new DoubleVecParam(InversionControlScalingFactorsEnum,control_scaling_factors,num_independent_objects));
     111       
     112                                /*cleanup*/
     113                                for(int i=0;i<num_independent_objects;i++){
     114                                        xDelete<char>(names[i]);
     115                                }
     116                                xDelete<char*>(names);
     117                                xDelete<int>(ind_enums);       
     118                                xDelete<int>(costfunc_enums);
     119                                break;
     120                        }
     121                        default:
     122                                _error_("not supported");
    59123                }
    60                 xDelete<char*>(cm_responses);
    61                 parameters->AddObject(new IntVecParam(InversionCostFunctionsEnum,costfunc_enums,num_costfunc));
    62                
    63                 xDelete<int>(control_enums);
    64                 xDelete<int>(costfunc_enums);
    65124
    66125                /*Inversion type specifics*/
    … …  
    107166                                parameters->AddObject(iomodel->CopyConstantObject("md.inversion.maxsteps",InversionMaxstepsEnum));
    108167                                parameters->AddObject(iomodel->CopyConstantObject("md.inversion.maxiter",InversionMaxiterEnum));
    109                                 iomodel->FetchData(&control_scaling_factors,NULL,NULL,"md.inversion.control_scaling_factors");
    110                                 parameters->AddObject(new DoubleVecParam(InversionControlScalingFactorsEnum,control_scaling_factors,num_controls));
    111                                 break;
     168                        break;
    112169                        default:
    113170                                _error_("not supported");
    114171                }
    115172
    116                 xDelete<int>(control_enums);
    117173                xDelete<int>(maxiter);
    118174                xDelete<IssmDouble>(control_scaling_factors);

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

    r22266 r22515  
    99#include "../ModelProcessorx.h"
    1010
     11
     12#if !defined(_HAVE_ADOLC_)
    1113void    UpdateElementsAndMaterialsControl(Elements* elements,Materials* materials, IoModel* iomodel){
    12 
     14       
    1315        /*Intermediary*/
    1416        bool       control_analysis;
    … …  
    2123        char     **cost_functions   = NULL;
    2224
     25        if(_HAVE_ADOLC_) return;
     26
    2327        /*Fetch parameters: */
    2428        iomodel->FindConstant(&control_analysis,"md.inversion.iscontrol");
    … …  
    2731        /*Now, return if no control*/
    2832        if(!control_analysis) return;
    29 
     33       
    3034        /*Process controls and convert from string to enums*/
    3135        iomodel->FindConstant(&controls,&num_controls,"md.inversion.control_parameters");
    … …  
    125129        xDelete<char*>(controls);
    126130}
     131#else
     132void UpdateElementsAndMaterialsControl(Elements* elements,Materials* materials, IoModel* iomodel){
     133
     134        /*Intermediaries*/
     135        int                             num_independent_objects,M,N;
     136        char**                  names                   = NULL;
     137        int*                            types                                                   = NULL;
     138        IssmDouble*             independent                                     = NULL;
     139        IssmDouble**    independents_min                        = NULL;
     140        IssmDouble**    independents_max                        = NULL;
     141
     142        /*Step1: create controls (independents)*/
     143        iomodel->FetchData(&num_independent_objects,"md.autodiff.num_independent_objects");
     144        _assert_(num_independent_objects>0);
     145        iomodel->FetchData(&names,&M,"md.autodiff.independent_object_names");
     146        _assert_(M==num_independent_objects);
     147        iomodel->FetchData(&types,NULL,NULL,"md.autodiff.independent_object_types");
     148
     149               
     150        /*TODO: fetch min and max*/
     151        independents_min = xNew<IssmDouble*>(num_independent_objects);
     152        independents_max = xNew<IssmDouble*>(num_independent_objects);
     153        for(int i=0;i<num_independent_objects;i++) independents_min[i]=NULL;
     154        for(int i=0;i<num_independent_objects;i++) independents_max[i]=NULL;
     155
     156        /*create independent objects, and at the same time, fetch the corresponding independent variables,
     157         *and declare them as such in ADOLC: */
     158        for(int i=0;i<num_independent_objects;i++){
     159
     160                if(types[i]==1){ /* vector:*/
     161
     162                        /*Get field name and input Enum from independent name*/
     163                        char* iofieldname  = NULL;
     164                        int   input_enum;
     165                        FieldAndEnumFromCode(&input_enum,&iofieldname,names[i]);
     166
     167                        /*Fetch required data*/
     168                        iomodel->FetchData(&independent,&M,&N,iofieldname);
     169                        iomodel->FetchData(&independents_min[i],&M,&N,"md.autodiff.independent_min_parameters");
     170                        iomodel->FetchData(&independents_max[i],&M,&N,"md.autodiff.independent_max_parameters");
     171                       
     172                        _assert_(independent);
     173
     174                        for(int j=0;j<elements->Size();j++){
     175                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
     176                                element->ControlInputCreate(independent,independents_min[i],independents_max[i],iomodel,M,N,input_enum,i+1);
     177                        }
     178                        xDelete<IssmDouble>(independent);
     179                }
     180                else{
     181                        _error_("Independent cannot be of size " << types[i]);
     182                }
     183        }
     184
     185        /*cleanup*/
     186        for(int i=0;i<num_independent_objects;i++){
     187                xDelete<char>(names[i]);
     188                xDelete<IssmDouble>(independents_min[i]);
     189                xDelete<IssmDouble>(independents_max[i]);
     190        }
     191        xDelete<IssmDouble*>(independents_min);
     192        xDelete<IssmDouble*>(independents_max);
     193        xDelete<char*>(names);
     194        xDelete<int>(types);
     195
     196        /*Step2: create cost functions (dependents)*/
     197
     198        return;
     199}
     200#endif

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

    r22323 r22515  
    1515        bool control_analysis;
    1616        bool dakota_analysis;
     17        bool adolc_analysis;
    1718        int nnat,dummy;
    1819        int* nature=NULL;
    … …  
    2223        iomodel->FindConstant(&dakota_analysis,"md.qmu.isdakota");
    2324        iomodel->FindConstant(&materials_type,"md.materials.type");
     25        iomodel->FindConstant(&adolc_analysis,"md.autodiff.isautodiff");
    2426
    2527        /*Did we already create the elements? : */
    … …  
    3032
    3133        /*Create elements*/
    32         if(control_analysis)iomodel->FetchData(2,"md.inversion.min_parameters","md.inversion.max_parameters");
     34        if(control_analysis && !adolc_analysis)iomodel->FetchData(2,"md.inversion.min_parameters","md.inversion.max_parameters");
     35       
    3336        switch(iomodel->meshelementtype){
    3437                case TriaEnum:
    … …  
    218221
    219222        /*Free data: */
    220         iomodel->DeleteData(7,"md.mesh.upperelements","md.mesh.lowerelements","md.material.rheology_B",
    221                                 "md.material.rheology_n","md.damage.D","md.inversion.min_parameters","md.inversion.max_parameters");
     223        iomodel->DeleteData(7,"md.mesh.upperelements","md.mesh.lowerelements","md.material.rheology_B","md.material.rheology_n","md.damage.D","md.inversion.min_parameters","md.inversion.max_parameters");
    222224
    223225        /*Add new constant material property to materials, at the end: */

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

    r22507 r22515  
    507507
    508508                                /*Intermediary*/
    509                                 int    numout;
    510                                 char **ncf_name_s             = NULL;
    511                                 char **ncf_definitionstring_s = NULL;
    512                                 char **cost_functions         = NULL;
    513                                 int    cost_function ,domaintype;
    514                                 int    num_cost_functions;
     509                                int          numout,numout2;
     510                                char       **ncf_name_s             = NULL;
     511                                char       **ncf_definitionstring_s = NULL;
     512                                char       **cost_functions         = NULL;
     513                                IssmDouble **cost_functions_weights = NULL;
     514                                int*         cost_functions_weights_M = NULL;
     515                                int*         cost_functions_weights_N = NULL;
     516                                int          cost_function,domaintype;
     517                                int          num_cost_functions;
    515518
    516519                                /*Process cost functions and convert from string to enums*/
    … …  
    523526                                }
    524527
    525                                 //iomodel->FetchData(1,"md.numberedcostfunction.cost_functions_coefficients");
    526 
    527528                                iomodel->FetchMultipleData(&ncf_name_s,&numout,"md.numberedcostfunction.name");
    528                                 iomodel->FetchMultipleData(&ncf_definitionstring_s,&numout,"md.numberedcostfunction.definitionstring");
    529                                 if(numout>1) _error_("not implemented yet, check code here");
     529                                iomodel->FetchMultipleData(&ncf_definitionstring_s,&numout2,"md.numberedcostfunction.definitionstring"); _assert_(numout2 == numout);
     530                                iomodel->FetchMultipleData(&cost_functions_weights,&cost_functions_weights_M,&cost_functions_weights_N,&numout2,"md.numberedcostfunction.cost_functions_coefficients");  _assert_(numout2 == numout);
     531                                if(numout!=1) _error_("not implemented yet, check code here");
    530532
    531533                                /*Fetch Observations */
    … …  
    543545                                                if(domaintype!=Domain2DverticalEnum) iomodel->FetchDataToInput(elements,"md.numberedcostfunction.vy_obs",InversionVyObsEnum);
    544546                                        }
    545 
    546547                                }
    547548
    548549                                for(j=0;j<numout;j++){
    549550
     551                                        /*Now, for this particular misfit object, make sure we plug into the elements: the observation, and the weights.*/
     552                                        for(int k=0;k<elements->Size();k++){
     553                                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(k));
     554                                                element->DatasetInputCreate(cost_functions_weights[i],cost_functions_weights_M[i],cost_functions_weights_N[i],cost_function_enums,num_cost_functions,iomodel,InversionCostFunctionsCoefficientsEnum);
     555                                        }
    550556                                        output_definitions->AddObject(new Numberedcostfunction(ncf_name_s[j],StringToEnumx(ncf_definitionstring_s[j]),num_cost_functions,cost_function_enums));
    551 
    552                                 }
    553 
     557                                }
    554558                               
    555559                                /*Free data: */
    556                                 //iomodel->DeleteData(1,"md.inversion.cost_functions_coefficients");
    557560                                iomodel->DeleteData(2,"md.numberedcostfunction.name","md.numberedcostfunction.definitionstring");
    558561                                xDelete<int>(cost_function_enums);
    … …  
    564567                                        xDelete<char>(ncf_name_s[j]);
    565568                                        xDelete<char>(ncf_definitionstring_s[j]);
     569                                        xDelete<IssmDouble>(cost_functions_weights[j]);
    566570                                }
    567571                                xDelete<char*>(ncf_name_s);
    568572                                xDelete<char*>(ncf_definitionstring_s);
     573                                xDelete<int>(cost_functions_weights_M);
     574                                xDelete<int>(cost_functions_weights_N);
     575                                xDelete<IssmDouble*>(cost_functions_weights);
    569576                        }
    570577                        /*}}}*/

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

    r22473 r22515  
    347347        #ifdef _HAVE_ADOLC_
    348348        if(VerboseMProcessor()) _printf0_("   starting autodiff parameters \n");
     349        printf("-------------- file: CreateParameters.cpp line: %i\n",__LINE__);
    349350        CreateParametersAutodiff(parameters,iomodel);
    350351        if(VerboseMProcessor()) _printf0_("   ending autodiff parameters \n");

  • issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp

    r22471 r22515  
    1111        /*output*/
    1212        char* fieldname = NULL;
    13         int   param_enum = -1;
     13        int   input_enum = -1;
    1414
    1515        if(strcmp(string_in,"Thickness")==0){
    1616                const char* field = "md.geometry.thickness";
    17                 param_enum        = ThicknessEnum;
     17                input_enum        = ThicknessEnum;
    1818                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    1919        }
    2020        else if(strcmp(string_in,"MaterialsRheologyB")==0){
    2121                const char* field = "md.materials.rheology_B";
    22                 param_enum        = MaterialsRheologyBEnum;
     22                input_enum        = MaterialsRheologyBEnum;
    2323                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    2424        }
    2525        else if(strcmp(string_in,"SmbMassBalance")==0){
    2626                const char* field = "md.smb.mass_balance";
    27                 param_enum        = SmbMassBalanceEnum;
     27                input_enum        = SmbMassBalanceEnum;
    2828                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    2929        }
    3030        else if(strcmp(string_in,"SmbAccumulation")==0){
    3131                const char* field = "md.smb.accumulation";
    32                 param_enum        = SmbAccumulationEnum;
     32                input_enum        = SmbAccumulationEnum;
    3333                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    3434        }
    3535        else if(strcmp(string_in,"SmbMelt")==0){
    3636                const char* field = "md.smb.melt";
    37                 param_enum        = SmbMeltEnum;
     37                input_enum        = SmbMeltEnum;
    3838                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    3939        }
    4040        else if(strcmp(string_in,"SmbRefreeze")==0){
    4141                const char* field = "md.smb.refreeze";
    42                 param_enum        = SmbRefreezeEnum;
     42                input_enum        = SmbRefreezeEnum;
    4343                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    4444        }
    4545        else if(strcmp(string_in,"SmbRunoff")==0){
    4646                const char* field = "md.smb.runoff";
    47                 param_enum        = SmbRunoffEnum;
     47                input_enum        = SmbRunoffEnum;
    4848                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    4949        }
    5050        else if(strcmp(string_in,"SmbEvaporation")==0){
    5151                const char* field = "md.smb.evaporation";
    52                 param_enum        = SmbEvaporationEnum;
     52                input_enum        = SmbEvaporationEnum;
    5353                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    5454        }
    5555        else if(strcmp(string_in,"SmbTa")==0){
    5656                const char* field = "md.smb.Ta";
    57                 param_enum        = SmbTaEnum;
     57                input_enum        = SmbTaEnum;
    5858                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    5959        }
    6060        else if(strcmp(string_in,"SmbV")==0){
    6161                const char* field = "md.smb.V";
    62                 param_enum        = SmbVEnum;
     62                input_enum        = SmbVEnum;
    6363                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    6464        }
    6565        else if(strcmp(string_in,"SmbDswrf")==0){
    6666                const char* field = "md.smb.dswrf";
    67                 param_enum        = SmbDswrfEnum;
     67                input_enum        = SmbDswrfEnum;
    6868                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    6969        }
    7070        else if(strcmp(string_in,"SmbDlwrf")==0){
    7171                const char* field = "md.smb.dlwrf";
    72                 param_enum        = SmbDlwrfEnum;
     72                input_enum        = SmbDlwrfEnum;
    7373                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    7474        }
    7575        else if(strcmp(string_in,"SmbP")==0){
    7676                const char* field = "md.smb.P";
    77                 param_enum        = SmbPEnum;
     77                input_enum        = SmbPEnum;
    7878                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    7979        }
    8080        else if(strcmp(string_in,"SmbEAir")==0){
    8181                const char* field = "md.smb.eAir";
    82                 param_enum        = SmbEAirEnum;
     82                input_enum        = SmbEAirEnum;
    8383                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    8484        }
    8585        else if(strcmp(string_in,"SmbPAir")==0){
    8686                const char* field = "md.smb.pAir";
    87                 param_enum        = SmbPAirEnum;
     87                input_enum        = SmbPAirEnum;
    8888                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    8989        }
    9090        else if(strcmp(string_in,"SmbVz")==0){
    9191                const char* field = "md.smb.Vz";
    92                 param_enum        = SmbVzEnum;
     92                input_enum        = SmbVzEnum;
    9393                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    9494        }
    9595        else if(strcmp(string_in,"SmbTz")==0){
    9696                const char* field = "md.smb.Tz";
    97                 param_enum        = SmbTzEnum;
     97                input_enum        = SmbTzEnum;
    9898                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    9999        }
    100100        else if(strcmp(string_in,"SmbC")==0){
    101101                const char* field = "md.smb.C";
    102                 param_enum        = SmbCEnum;
     102                input_enum        = SmbCEnum;
    103103                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    104104        }
    105105        else if(strcmp(string_in,"BasalforcingsFloatingiceMeltingRate")==0){
    106106                const char* field = "md.basalforcings.floatingice_melting_rate";
    107                 param_enum        = BasalforcingsFloatingiceMeltingRateEnum;
    108                 fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    109         }
    110         else if(strcmp(string_in,"FrictionCoefficient")==0){
     107                input_enum        = BasalforcingsFloatingiceMeltingRateEnum;
     108                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
     109        }
     110        else if(strcmp(string_in,"FrictionCoefficient")==0 || strcmp(string_in,"md.friction.coefficient")==0){
    111111                const char* field = "md.friction.coefficient";
    112                 param_enum        = FrictionCoefficientEnum;
     112                input_enum        = FrictionCoefficientEnum;
    113113                fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    114114        }
    115115        else if(strcmp(string_in,"Vx")==0){
    116116                 const char* field = "md.initialization.vx";
    117                  param_enum        = VxEnum;
     117                 input_enum        = VxEnum;
    118118                 fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    119119         }
    120120        else if(strcmp(string_in,"Vy")==0){
    121121                 const char* field = "md.initialization.vy";
    122                  param_enum        = VyEnum;
     122                 input_enum        = VyEnum;
    123123                 fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    124124         }
    125125        else if(strcmp(string_in,"BalancethicknessThickeningRate")==0){
    126126                 const char* field = "md.balancethickness.thickening_rate";
    127                  param_enum        = BalancethicknessThickeningRateEnum;
     127                 input_enum        = BalancethicknessThickeningRateEnum;
    128128                 fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    129129        }
    130130        else if(strcmp(string_in,"BalancethicknessSpcthickness")==0){
    131131                 const char* field = "md.balancethickness.spcthickness";
    132                  param_enum        = BalancethicknessSpcthicknessEnum;
     132                 input_enum        = BalancethicknessSpcthicknessEnum;
    133133                 fieldname=xNew<char>((strlen(field)+1)); xMemCpy<char>(fieldname,field,(strlen(field)+1));
    134134        }
    … …  
    138138
    139139        /*Assign output pointers*/
    140         *out_enum = param_enum;
     140        *out_enum = input_enum;
    141141        *pfield   = fieldname;
    142142        return;

  • issm/trunk-jpl/src/m/classes/adm1qn3inversion.m

    r22443 r22515  
    77        properties (SetAccess=public)
    88                iscontrol                   = 0
    9                 control_parameters          = NaN
    10                 control_scaling_factors     = NaN
    119                maxsteps                    = 0
    1210                maxiter                     = 0
    1311                dxmin                       = 0
    1412                gttol                       = 0
    15                 cost_functions              = NaN
    16                 cost_functions_coefficients = NaN
    17                 min_parameters              = NaN
    18                 max_parameters              = NaN
    19                 vx_obs                      = NaN
    20                 vy_obs                      = NaN
    21                 vz_obs                      = NaN
    22                 vel_obs                     = NaN
    23                 thickness_obs               = NaN
    24                 surface_obs                 = NaN
    2513
    2614        end
    … …  
    3725                end % }}}
    3826                function self = extrude(self,md) % {{{
    39                         self.vx_obs=project3d(md,'vector',self.vx_obs,'type','node');
    40                         self.vy_obs=project3d(md,'vector',self.vy_obs,'type','node');
    41                         self.vel_obs=project3d(md,'vector',self.vel_obs,'type','node');
    42                         self.thickness_obs=project3d(md,'vector',self.thickness_obs,'type','node');
    43                         if numel(self.cost_functions_coefficients)>1,self.cost_functions_coefficients=project3d(md,'vector',self.cost_functions_coefficients,'type','node');end;
    44                         if numel(self.min_parameters)>1,self.min_parameters=project3d(md,'vector',self.min_parameters,'type','node');end;
    45                         if numel(self.max_parameters)>1,self.max_parameters=project3d(md,'vector',self.max_parameters,'type','node');end;
    4627                end % }}}
    4728                function self = setdefaultparameters(self) % {{{
    … …  
    5031                        %parameter to be inferred by control methods (only
    5132                        %drag and B are supported yet)
    52                         self.control_parameters={'FrictionCoefficient'};
    53 
    54                         %Scaling factor for each control
    55                         self.control_scaling_factors=1;
     33                        %self.control_parameters={'FrictionCoefficient'};
    5634
    5735                        %number of iterations
    5836                        self.maxsteps=20;
    5937                        self.maxiter=40;
    60 
    61                         %several responses can be used:
    62                         self.cost_functions=101;
    6338
    6439                        %m1qn3 parameters
    … …  
    7550                                md = checkmessage(md,['M1QN3 has not been installed, ISSM needs to be reconfigured and recompiled with M1QN3']);
    7651                        end
    77                         num_controls=numel(md.inversion.control_parameters);
    78                         num_costfunc=size(md.inversion.cost_functions,2);
    7952
    8053                        md = checkfield(md,'fieldname','inversion.iscontrol','values',[0 1]);
    81                         md = checkfield(md,'fieldname','inversion.control_parameters','cell',1,'values',supportedcontrols());
    82                         md = checkfield(md,'fieldname','inversion.control_scaling_factors','size',[1 num_controls],'>',0,'NaN',1,'Inf',1);
    8354                        md = checkfield(md,'fieldname','inversion.maxsteps','numel',1,'>=',0);
    8455                        md = checkfield(md,'fieldname','inversion.maxiter','numel',1,'>=',0);
    8556                        md = checkfield(md,'fieldname','inversion.dxmin','numel',1,'>',0);
    8657                        md = checkfield(md,'fieldname','inversion.gttol','numel',1,'>',0);
    87                         md = checkfield(md,'fieldname','inversion.cost_functions','size',[1 num_costfunc],'values',supportedcostfunctions());
    88                         md = checkfield(md,'fieldname','inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices num_costfunc],'>=',0);
    89                         md = checkfield(md,'fieldname','inversion.min_parameters','size',[md.mesh.numberofvertices num_controls]);
    90                         md = checkfield(md,'fieldname','inversion.max_parameters','size',[md.mesh.numberofvertices num_controls]);
    9158
    92                         if strcmp(solution,'BalancethicknessSolution')
    93                                 md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1);
    94                                 md = checkfield(md,'fieldname','inversion.surface_obs','size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1);
    95                         elseif strcmp(solution,'BalancethicknessSoftSolution')
    96                                 md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1);
    97                         else
    98                                 md = checkfield(md,'fieldname','inversion.vx_obs','size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1);
    99                                 if ~strcmp(domaintype(md.mesh),'2Dvertical'),
    100                                         md = checkfield(md,'fieldname','inversion.vy_obs','size',[md.mesh.numberofvertices 1],'NaN',1,'Inf',1);
    101                                 end
    102                         end
    103                 end % }}}
     59        end % }}}
    10460                function disp(self) % {{{
    10561                        disp(sprintf('   adm1qn3inversion parameters:'));
    10662                        fielddisplay(self,'iscontrol','is inversion activated?');
    107                         fielddisplay(self,'control_parameters','ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''}');
    108                         fielddisplay(self,'control_scaling_factors','order of magnitude of each control (useful for multi-parameter optimization)');
    10963                        fielddisplay(self,'maxsteps','maximum number of iterations (gradient computation)');
    11064                        fielddisplay(self,'maxiter','maximum number of Function evaluation (forward run)');
    11165                        fielddisplay(self,'dxmin','convergence criterion: two points less than dxmin from eachother (sup-norm) are considered identical');
    11266                        fielddisplay(self,'gttol','convergence criterion: ||g(X)||/||g(X0)|| (g(X0): gradient at initial guess X0)');
    113                         fielddisplay(self,'cost_functions','indicate the type of response for each optimization step');
    114                         fielddisplay(self,'cost_functions_coefficients','cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter');
    115                         fielddisplay(self,'min_parameters','absolute minimum acceptable value of the inversed parameter on each vertex');
    116                         fielddisplay(self,'max_parameters','absolute maximum acceptable value of the inversed parameter on each vertex');
    117                         fielddisplay(self,'vx_obs','observed velocity x component [m/yr]');
    118                         fielddisplay(self,'vy_obs','observed velocity y component [m/yr]');
    119                         fielddisplay(self,'vel_obs','observed velocity magnitude [m/yr]');
    120                         fielddisplay(self,'thickness_obs','observed thickness [m]');
    121                         fielddisplay(self,'surface_obs','observed surface elevation [m]');
    122                         disp('Available cost functions:');
    123                         disp('   101: SurfaceAbsVelMisfit');
    124                         disp('   102: SurfaceRelVelMisfit');
    125                         disp('   103: SurfaceLogVelMisfit');
    126                         disp('   104: SurfaceLogVxVyMisfit');
    127                         disp('   105: SurfaceAverageVelMisfit');
    128                         disp('   201: ThicknessAbsMisfit');
    129                         disp('   501: DragCoefficientAbsGradient');
    130                         disp('   502: RheologyBbarAbsGradient');
    131                         disp('   503: ThicknessAbsGradient');
    13267                end % }}}
    13368                function marshall(self,prefix,md,fid) % {{{
    … …  
    13873                        WriteData(fid,prefix,'name','md.inversion.type','data',4,'format','Integer');
    13974                        if ~self.iscontrol, return; end
    140                         WriteData(fid,prefix,'object',self,'class','inversion','fieldname','control_scaling_factors','format','DoubleMat','mattype',3);
    14175                        WriteData(fid,prefix,'object',self,'class','inversion','fieldname','maxsteps','format','Integer');
    14276                        WriteData(fid,prefix,'object',self,'class','inversion','fieldname','maxiter','format','Integer');
    14377                        WriteData(fid,prefix,'object',self,'class','inversion','fieldname','dxmin','format','Double');
    14478                        WriteData(fid,prefix,'object',self,'class','inversion','fieldname','gttol','format','Double');
    145                         WriteData(fid,prefix,'object',self,'class','inversion','fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1);
    146                         WriteData(fid,prefix,'object',self,'class','inversion','fieldname','min_parameters','format','DoubleMat','mattype',3);
    147                         WriteData(fid,prefix,'object',self,'class','inversion','fieldname','max_parameters','format','DoubleMat','mattype',3);
    148                         WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vx_obs','format','DoubleMat','mattype',1,'scale',1./yts);
    149                         WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vy_obs','format','DoubleMat','mattype',1,'scale',1./yts);
    150                         WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vz_obs','format','DoubleMat','mattype',1,'scale',1./yts);
    151                         WriteData(fid,prefix,'object',self,'class','inversion','fieldname','vel_obs','format','DoubleMat','mattype',1,'scale',1./yts);
    152                         if(numel(self.thickness_obs)==md.mesh.numberofelements),
    153                                 mattype=2;
    154                         else
    155                                 mattype=1;
    156                         end
    157                         WriteData(fid,prefix,'object',self,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',mattype);
    158                         WriteData(fid,prefix,'object',self,'class','inversion','fieldname','surface_obs','format','DoubleMat','mattype',mattype);
    15979
    160                         %process control parameters
    161                         num_control_parameters=numel(self.control_parameters);
    162                         WriteData(fid,prefix,'object',self,'fieldname','control_parameters','format','StringArray');
    163                         WriteData(fid,prefix,'data',num_control_parameters,'name','md.inversion.num_control_parameters','format','Integer');
    164 
    165                         %process cost functions
    166                         num_cost_functions=size(self.cost_functions,2);
    167                         data=marshallcostfunctions(self.cost_functions);
    168                         WriteData(fid,prefix,'data',data,'name','md.inversion.cost_functions','format','StringArray');
    169                         WriteData(fid,prefix,'data',num_cost_functions,'name','md.inversion.num_cost_functions','format','Integer');
    17080                end % }}}
    17181                function savemodeljs(self,fid,modelname) % {{{
    17282               
    17383                        writejsdouble(fid,[modelname '.inversion.iscontrol'],self.iscontrol);
    174                         writejscellstring(fid,[modelname '.inversion.control_parameters'],self.control_parameters);
    175                         writejsdouble(fid,[modelname '.inversion.control_scaling_factors'],self.control_scaling_factors);
    17684                        writejsdouble(fid,[modelname '.inversion.maxsteps'],self.maxsteps);
    17785                        writejsdouble(fid,[modelname '.inversion.maxiter'],self.maxiter);
    17886                        writejsdouble(fid,[modelname '.inversion.dxmin'],self.dxmin);
    17987                        writejsdouble(fid,[modelname '.inversion.gttol'],self.gttol);
    180                         writejs2Darray(fid,[modelname '.inversion.cost_functions'],self.cost_functions);
    181                         writejs2Darray(fid,[modelname '.inversion.cost_functions_coefficients'],self.cost_functions_coefficients);
    182                         writejs1Darray(fid,[modelname '.inversion.min_parameters'],self.min_parameters);
    183                         writejs1Darray(fid,[modelname '.inversion.max_parameters'],self.max_parameters);
    184                         writejs1Darray(fid,[modelname '.inversion.vx_obs'],self.vx_obs);
    185                         writejs1Darray(fid,[modelname '.inversion.vy_obs'],self.vy_obs);
    186                         writejs1Darray(fid,[modelname '.inversion.vz_obs'],self.vz_obs);
    187                         writejs1Darray(fid,[modelname '.inversion.vel_obs'],self.vel_obs);
    188                         writejs1Darray(fid,[modelname '.inversion.thickness_obs'],self.thickness_obs);
    189                         writejs1Darray(fid,[modelname '.inversion.surface_obs'],self.surface_obs);
    19088
    19189                end % }}}

  • issm/trunk-jpl/src/m/classes/autodiff.m

    r20979 r22515  
    130130                                        names{i}=indep.name;
    131131                                        types(i)=indep.typetoscalar();
     132                                        min_parameters(:,i)=indep.min_parameters;
     133                                        max_parameters(:,i)=indep.max_parameters;
     134                                        scaling_factors(i)=indep.control_scaling_factor;
    132135                                end
    133136                                WriteData(fid,prefix,'data',names,'name','md.autodiff.independent_object_names','format','StringArray');
    134137                                WriteData(fid,prefix,'data',types,'name','md.autodiff.independent_object_types','format','IntMat','mattype',3);
     138                                WriteData(fid,prefix,'data',min_parameters,'name','md.autodiff.independent_min_parameters','format','DoubleMat','mattype',3);
     139                 WriteData(fid,prefix,'data',max_parameters,'name','md.autodiff.independent_max_parameters','format','DoubleMat','mattype',3);
     140                 WriteData(fid,prefix,'data',scaling_factors,'name','md.autodiff.independent_scaling_factors','format','DoubleMat','mattype',3);
     141
     142
    135143                        end
    136144                        %}}}

  • issm/trunk-jpl/src/m/classes/independent.m

    r21049 r22515  
    1111                fov_forward_indices  = [];
    1212                nods                 = 0;
     13                min_parameters                  = NaN;
     14                max_parameters                  = NaN;
     15                control_scaling_factor     = NaN
    1316        end
    1417        methods
    … …  
    4548                                end
    4649                                md = checkfield(md,'fieldname',['autodiff.independents{' num2str(i) '}.fov_forward_indices'],'>=',1,'<=',self.nods,'size',[NaN 1]);
     50                                md = checkfield(md,'fieldname',['autodiff.independents{' num2str(i) '}.min_parameters'],'size',[md.mesh.numberofvertices 1]);
     51                                md = checkfield(md,'fieldname',['autodiff.independents{' num2str(i) '}.max_parameters'],'size',[md.mesh.numberofvertices 1]);
     52                                md = checkfield(md,'fieldname',['autodiff.independents{' num2str(i) '}.control_scaling_factors'],'size',[1 1],'>',0,'NaN',1,'Inf',1);
     53
    4754                        end
    4855
    … …  
    5461                        fielddisplay(self,'type','type of variable (''vertex'' or ''scalar'')');
    5562                        fielddisplay(self,'nods','size of dependent variables');
     63                        fielddisplay(self,'min_parameters','absolute minimum acceptable value of the inversed parameter on each vertex');
     64                        fielddisplay(self,'max_parameters','absolute maximum acceptable value of the inversed parameter on each vertex');
     65                        fielddisplay(self,'control_scaling_factor','order of magnitude of each control (useful for multi-parameter optimization)');
    5666                        if ~isnan(self.fos_forward_index),
    5767                                fielddisplay(self,'fos_forward_index','index for fos_foward driver of ADOLC');
Note: See TracChangeset for help on using the changeset viewer.