Changeset 4967


Ignore:
Timestamp:
08/04/10 11:01:57 (15 years ago)
Author:
Mathieu Morlighem
Message:

Update control method: now use module InputControlUpdate to update/constrain and save the parameter

Location:
issm/trunk/src
Files:
6 added
23 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Container/Inputs.cpp

    r4943 r4967  
    426426}
    427427/*}}}*/
     428/*FUNCTION Inputs::AXPY{{{1*/
     429void  Inputs::AXPY(int YEnum, double scalar, int XEnum){
     430           
     431        Input* xinput=NULL;
     432        Input* yinput=NULL;
     433
     434        /*Find x and y inputs: */
     435        xinput=(Input*)this->GetInput(XEnum);
     436        yinput=(Input*)this->GetInput(YEnum);
     437
     438        /*some checks: */
     439        if(!xinput) ISSMERROR(" input %s could not be found!",EnumAsString(XEnum));
     440        if(!yinput) ISSMERROR(" input %s could not be found!",EnumAsString(YEnum));
     441
     442        /*Apply AXPY: */
     443        yinput->AXPY(xinput,scalar);
     444}
     445/*}}}*/

  • issm/trunk/src/c/Container/Inputs.h

    r4943 r4967  
    3434                Input*  GetInput(int enum_name);
    3535                Inputs* SpawnTriaInputs(int* indices);
     36                void    AXPY(int YEnum, double scalar, int XEnum);
    3637               
    3738                void GetParameterAverage(double* pvalue, int enum_type);

  • issm/trunk/src/c/Makefile.am

    r4939 r4967  
    435435                                        ./modules/InputControlConstrainx/InputControlConstrainx.h\
    436436                                        ./modules/InputControlConstrainx/InputControlConstrainx.cpp\
     437                                        ./modules/InputControlUpdatex/InputControlUpdatex.h\
     438                                        ./modules/InputControlUpdatex/InputControlUpdatex.cpp\
    437439                                        ./modules/SurfaceAreax/SurfaceAreax.h\
    438440                                        ./modules/SurfaceAreax/SurfaceAreax.cpp\
    … …  
    978980                                        ./modules/InputControlConstrainx/InputControlConstrainx.h\
    979981                                        ./modules/InputControlConstrainx/InputControlConstrainx.cpp\
     982                                        ./modules/InputControlUpdatex/InputControlUpdatex.h\
     983                                        ./modules/InputControlUpdatex/InputControlUpdatex.cpp\
    980984                                        ./modules/SurfaceAreax/SurfaceAreax.h\
    981985                                        ./modules/SurfaceAreax/SurfaceAreax.cpp\

  • issm/trunk/src/c/modules/Gradjx/Gradjx.cpp

    r4592 r4967  
    1515        int  dim;
    1616        int  numberofvertices;
    17         bool extrude;
    1817        Vec  gradient = NULL;
    1918       
    2019        /*retrieve some parameters: */
    2120        parameters->FindParam(&dim,DimEnum);
    22         parameters->FindParam(&extrude,ExtrudeParamEnum);
    2321        numberofvertices=vertices->NumberOfVertices();
    2422
    … …  
    3634        VecAssemblyEnd(gradient);
    3735
    38         /*Extrude if needed: */
    39         if(dim==3 && extrude) VecExtrudex(gradient, elements,nodes, vertices, loads, materials, parameters,0);
    40 
    4136        /*Assign output pointers: */
    4237        *pgradient=gradient;

  • issm/trunk/src/c/modules/InputDuplicatex/InputDuplicatex.cpp

    r4572 r4967  
    1212       
    1313        /*Go through elemnets, and ask to reinitialie the input: */
    14 
    1514        int      i;
    1615        for(i=0;i<elements->Size();i++){
    … …  
    1817                element->InputDuplicate(original_enum,new_enum);
    1918        }
     19        for(i=0;i<materials->Size();i++){
     20                Material* material=(Material*)materials->GetObjectByOffset(i);
     21                //material->InputDuplicate(original_enum,new_enum);
     22        }
    2023
    2124}

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

    r4919 r4967  
    4040        IoModelFetchData(&iomodel->rheology_B,NULL,NULL,iomodel_handle,"rheology_B");
    4141        IoModelFetchData(&iomodel->rheology_n,NULL,NULL,iomodel_handle,"rheology_n");
     42        if(iomodel->control_analysis && strcmp(iomodel->control_type,"rheology_B")==0){
     43                IoModelFetchData(&iomodel->control_parameter,NULL,NULL,iomodel_handle,iomodel->control_type); //copy the control parameter in iomodel
     44        }
    4245       
    4346        /*Create elements and materials: */
    … …  
    6164        xfree((void**)&iomodel->rheology_B);
    6265        xfree((void**)&iomodel->rheology_n);
     66        xfree((void**)&iomodel->control_parameter);
    6367
    6468        /*Add new constrant material property tgo materials, at the end: */

  • issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp

    r4578 r4967  
    4747                IoModelFetchData(&iomodel->vy_obs,NULL,NULL,iomodel_handle,"vy_obs");
    4848                IoModelFetchData(&iomodel->weights,NULL,NULL,iomodel_handle,"weights");
    49                 IoModelFetchData(&iomodel->control_parameter,NULL,NULL,iomodel_handle,iomodel->control_type); //copy the control parameter in iomodel
     49                if(strcmp(iomodel->control_type,"drag_coefficient")==0){
     50                        IoModelFetchData(&iomodel->control_parameter,NULL,NULL,iomodel_handle,iomodel->control_type); //copy the control parameter in iomodel
     51                }
    5052        }
    5153

  • issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/UpdateElementsDiagnosticStokes.cpp

    r4578 r4967  
    4848                IoModelFetchData(&iomodel->vy_obs,NULL,NULL,iomodel_handle,"vy_obs");
    4949                IoModelFetchData(&iomodel->weights,NULL,NULL,iomodel_handle,"weights");
    50                 IoModelFetchData(&iomodel->control_parameter,NULL,NULL,iomodel_handle,iomodel->control_type); //copy the control parameter in iomodel
     50                if(strcmp(iomodel->control_type,"drag_coefficient")==0){
     51                        IoModelFetchData(&iomodel->control_parameter,NULL,NULL,iomodel_handle,iomodel->control_type); //copy the control parameter in iomodel
     52                }
    5153        }
    5254

  • issm/trunk/src/c/modules/modules.h

    r4939 r4967  
    2828#include "./InputAXPYx/InputAXPYx.h"
    2929#include "./InputControlConstrainx/InputControlConstrainx.h"
     30#include "./InputControlUpdatex/InputControlUpdatex.h"
    3031#include "./InputConvergencex/InputConvergencex.h"
    3132#include "./InputDuplicatex/InputDuplicatex.h"

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

    r4942 r4967  
    7272                virtual void   GetVectorFromInputs(Vec vector,int NameEnum)=0;
    7373                virtual void   InputAXPY(int YEnum, double scalar, int XEnum)=0;
     74                virtual void   InputControlUpdate(double scalar,bool save_parameter)=0;
    7475                virtual void   InputControlConstrain(int control_type,double cm_min, double cm_max)=0;
    7576                virtual bool   InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums)=0;

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

    r4945 r4967  
    10261026}
    10271027/*}}}*/
     1028/*FUNCTION Penta::InputControlUpdate{{{1*/
     1029void  Penta::InputControlUpdate(double scalar,bool save_parameter){
     1030
     1031        /*Intermediary*/
     1032        Input* input=NULL;
     1033        double cm_min,cm_max;
     1034        int    control_type;
     1035
     1036        /*retrieve some parameters: */
     1037        this->parameters->FindParam(&cm_min,CmMinEnum);
     1038        this->parameters->FindParam(&cm_max,CmMaxEnum);
     1039        this->parameters->FindParam(&control_type,ControlTypeEnum);
     1040
     1041
     1042        /*Rheology*/
     1043        if(control_type==RheologyB2dEnum){
     1044
     1045                /*First, get revert to previous parameter value (kept in ControlParameter input)*/
     1046                matice->inputs->DuplicateInput(ControlParameterEnum,RheologyBEnum);
     1047
     1048                /*Now, update using the gradient new = old + scalar * gradient*/
     1049                //matice->inputs->AXPY(RheologyB2dEnum,scalar,GradientEnum);
     1050                // For now: Gradient is in element (TO BE CHANGED) and parameter in matice
     1051                Input* input_B   =(Input*)matice->inputs->GetInput(RheologyBEnum); ISSMASSERT(input_B);
     1052                Input* input_Grad=(Input*)this->inputs->GetInput(GradientEnum); ISSMASSERT(input_Grad);
     1053                input_B->AXPY(input_Grad,scalar);
     1054
     1055                /*Constrain constrain input*/
     1056                input=(Input*)matice->inputs->GetInput(RheologyBEnum); ISSMASSERT(input);
     1057                input->Constrain(cm_min,cm_max);
     1058
     1059                /*Finally: save input if requested*/
     1060                if (save_parameter) matice->inputs->DuplicateInput(RheologyBEnum,ControlParameterEnum);
     1061
     1062        }
     1063        else if(control_type==DragCoefficientEnum){
     1064
     1065                /*First, get revert to previous parameter value (kept in ControlParameter input)*/
     1066                this->inputs->DuplicateInput(ControlParameterEnum,DragCoefficientEnum);
     1067
     1068                /*Now, update using the gradient new = old + scalar * gradient*/
     1069                this->inputs->AXPY(DragCoefficientEnum,scalar,GradientEnum);
     1070
     1071                /*Constrain input*/
     1072                input=(Input*)this->inputs->GetInput(DragCoefficientEnum); ISSMASSERT(input);
     1073                input->Constrain(cm_min,cm_max);
     1074
     1075                /*Finally: save input if requested*/
     1076                if (save_parameter) inputs->DuplicateInput(DragCoefficientEnum,ControlParameterEnum);
     1077
     1078        }
     1079        else{
     1080                ISSMERROR("control type %s not implemented yet",EnumAsString(control_type));
     1081        }
     1082}
     1083/*}}}*/
    10281084/*FUNCTION Penta::InputDepthAverageAtBase{{{1*/
    10291085void  Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum){
    … …  
    11561212
    11571213        /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */
    1158         for (i=0;i<this->inputs->Size();i++){
    1159                 input=(Input*)this->inputs->GetObjectByOffset(i);
    1160                 if (input->EnumType()==enum_type){
    1161                         found=true; break;
    1162                 }
    1163         }
    1164         if (!found) ISSMERROR("Input %s not found in penta->inputs",EnumAsString(enum_type));
     1214        if (enum_type==RheologyB2dEnum){
     1215                input=this->matice->inputs->GetInput(enum_type);
     1216        }
     1217        else{
     1218                input=this->inputs->GetInput(enum_type);
     1219        }
     1220        if (!input) ISSMERROR("Input %s not found in tria->inputs",EnumAsString(enum_type));
    11651221
    11661222        /*If we don't find it, no big deal, just don't do the transfer. Otherwise, build a new Result
    … …  
    17481804                this->inputs->AddInput(new IntInput(DragTypeEnum,iomodel->drag_type));
    17491805
    1750         }
    1751         if (iomodel->rheology_B) {
    1752                 for(i=0;i<6;i++)nodeinputs[i]=iomodel->rheology_B[penta_vertex_ids[i]-1];
    1753                 this->inputs->AddInput(new PentaVertexInput(RheologyBEnum,nodeinputs));
    17541806        }
    17551807        if (iomodel->control_parameter) {
    … …  
    53925444                                name==VzEnum ||
    53935445                                name==TemperatureEnum ||
    5394                                 name==RheologyBEnum ||
    5395                                 name==RheologyNEnum ||
     5446                                name==ControlParameterEnum ||
    53965447                                name==FitEnum ||
    53975448                                name==DragCoefficientEnum ||

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

    r4942 r4967  
    8787                void   GradjDrag(Vec gradient);
    8888                void   InputAXPY(int YEnum, double scalar, int XEnum);
     89                void   InputControlUpdate(double scalar,bool save_parameter);
    8990                void   InputControlConstrain(int control_type,double cm_min, double cm_max);
    9091                bool   InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums);

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

    r4943 r4967  
    416416
    417417                        /*update input*/
    418                         this->inputs->AddInput(new TriaVertexInput(name,values));
     418                        if (name==RheologyB2dEnum || name==RheologyBEnum){
     419                                matice->inputs->AddInput(new TriaVertexInput(name,values));
     420                        }
     421                        else{
     422                                this->inputs->AddInput(new TriaVertexInput(name,values));
     423                        }
    419424                        return;
    420425
    … …  
    628633                }
    629634                else if (control_type==RheologyB2dEnum){
    630                         inputs->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], &gauss_l1l2l3[0],RheologyB2dEnum);
     635                        matice->inputs->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], &gauss_l1l2l3[0],RheologyB2dEnum);
    631636                        Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight;
    632637                }
    … …  
    965970        adjointx_input=inputs->GetInput(AdjointxEnum);
    966971        adjointy_input=inputs->GetInput(AdjointyEnum);
    967         rheologyb_input=inputs->GetInput(RheologyB2dEnum);
     972        rheologyb_input=matice->inputs->GetInput(RheologyB2dEnum);
    968973
    969974        /* Start  looping on the number of gaussian points: */
    … …  
    12291234        yinput->AXPY(xinput,scalar);
    12301235
    1231         /*Move input to Material if required (needed if control method) TO BE IMPROVED*/
    1232         if (YEnum==RheologyBEnum || YEnum==RheologyB2dEnum){
    1233                 this->matice->inputs->AddInput((Input*)yinput->copy());
     1236}
     1237/*}}}*/
     1238/*FUNCTION Tria::InputControlUpdate{{{1*/
     1239void  Tria::InputControlUpdate(double scalar,bool save_parameter){
     1240
     1241        /*Intermediary*/
     1242        Input* input=NULL;
     1243        double cm_min,cm_max;
     1244        int    control_type;
     1245
     1246        /*retrieve some parameters: */
     1247        this->parameters->FindParam(&cm_min,CmMinEnum);
     1248        this->parameters->FindParam(&cm_max,CmMaxEnum);
     1249        this->parameters->FindParam(&control_type,ControlTypeEnum);
     1250
     1251        /*Rheology*/
     1252        if(control_type==RheologyB2dEnum){
     1253
     1254                /*First, get revert to previous parameter value (kept in ControlParameter input)*/
     1255                matice->inputs->DuplicateInput(ControlParameterEnum,RheologyB2dEnum);
     1256
     1257                /*Now, update using the gradient new = old + scalar * gradient*/
     1258                //matice->inputs->AXPY(RheologyB2dEnum,scalar,GradientEnum);
     1259                // For now: Gradient is in element (TO BE CHANGED) and parameter in matice
     1260                Input* input_B   =(Input*)matice->inputs->GetInput(RheologyB2dEnum); ISSMASSERT(input_B);
     1261                Input* input_Grad=(Input*)this->inputs->GetInput(GradientEnum); ISSMASSERT(input_Grad);
     1262                input_B->AXPY(input_Grad,scalar);
     1263
     1264                /*Constrain constrain input*/
     1265                input=(Input*)matice->inputs->GetInput(RheologyB2dEnum); ISSMASSERT(input);
     1266                input->Constrain(cm_min,cm_max);
     1267
     1268                /*Finally: save input if requested*/
     1269                if (save_parameter) matice->inputs->DuplicateInput(RheologyB2dEnum,ControlParameterEnum);
     1270
     1271        }
     1272        else if(control_type==DragCoefficientEnum){
     1273
     1274                /*First, get revert to previous parameter value (kept in ControlParameter input)*/
     1275                this->inputs->DuplicateInput(ControlParameterEnum,DragCoefficientEnum);
     1276
     1277                /*Now, update using the gradient new = old + scalar * gradient*/
     1278                this->inputs->AXPY(DragCoefficientEnum,scalar,GradientEnum);
     1279
     1280                /*Constrain input*/
     1281                input=(Input*)this->inputs->GetInput(DragCoefficientEnum); ISSMASSERT(input);
     1282                input->Constrain(cm_min,cm_max);
     1283
     1284                /*Finally: save input if requested*/
     1285                if (save_parameter) inputs->DuplicateInput(DragCoefficientEnum,ControlParameterEnum);
     1286
     1287        }
     1288        else{
     1289                ISSMERROR("control type %s not implemented yet",EnumAsString(control_type));
    12341290        }
    12351291
    … …  
    13191375
    13201376        /*Call inputs method*/
    1321         inputs->DuplicateInput(original_enum,new_enum);
     1377        if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum);
    13221378
    13231379}
    … …  
    13441400
    13451401        /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */
    1346         for (i=0;i<this->inputs->Size();i++){
    1347                 input=(Input*)this->inputs->GetObjectByOffset(i);
    1348                 if (input->EnumType()==enum_type){
    1349                         found=true;
    1350                         break;
    1351                 }
    1352         }
    1353         if (!found) ISSMERROR("Input %s not found in tria->inputs",EnumAsString(enum_type));
     1402        if (enum_type==RheologyB2dEnum){
     1403                input=this->matice->inputs->GetInput(enum_type);
     1404        }
     1405        else{
     1406                input=this->inputs->GetInput(enum_type);
     1407        }
     1408        if (!input) ISSMERROR("Input %s not found in tria->inputs",EnumAsString(enum_type));
    13541409
    13551410        /*If we don't find it, no big deal, just don't do the transfer. Otherwise, build a new Result
    … …  
    21222177                if (iomodel->drag_q) this->inputs->AddInput(new DoubleInput(DragQEnum,iomodel->drag_q[index]));
    21232178                this->inputs->AddInput(new IntInput(DragTypeEnum,iomodel->drag_type));
    2124         }
    2125         if (iomodel->rheology_B) {
    2126                 for(i=0;i<3;i++)nodeinputs[i]=iomodel->rheology_B[tria_vertex_ids[i]-1];
    2127                 this->inputs->AddInput(new TriaVertexInput(RheologyB2dEnum,nodeinputs));
    21282179        }
    21292180        if (iomodel->control_parameter) {
    … …  
    60446095                                name==MeltingRateEnum ||
    60456096                                name==AccumulationRateEnum ||
     6097                                name==ControlParameterEnum ||
    60466098                                name==VxEnum ||
    60476099                                name==VyEnum ||
    6048                                 name==RheologyBEnum ||
    6049                                 name==RheologyB2dEnum ||
    6050                                 name==RheologyNEnum ||
    60516100                                name==FitEnum ||
    60526101                                name==DragCoefficientEnum ||

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

    r4942 r4967  
    8383                void   GradjDrag(Vec gradient);
    8484                void   InputAXPY(int YEnum, double scalar, int XEnum);
     85                void   InputControlUpdate(double scalar,bool save_parameter);
    8586                void   InputControlConstrain(int control_type,double cm_min, double cm_max);
    8687                bool   InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums);

  • issm/trunk/src/c/objects/Materials/Material.h

    r4248 r4967  
    1818        public:
    1919                virtual       ~Material(){};
     20
     21                /*Numerics*/
     22                virtual void   InputDuplicate(int original_enum,int new_enum)=0;
     23
    2024};
    2125#endif

  • issm/trunk/src/c/objects/Materials/Matice.cpp

    r4931 r4967  
    5555                        this->inputs->AddInput(new TriaVertexInput(RheologyNEnum,nodeinputs));
    5656                }
     57
     58                /*Get control_parameter*/
     59                if (iomodel->control_parameter) {
     60                        for(i=0;i<num_vertices;i++)nodeinputs[i]=iomodel->control_parameter[int(iomodel->elements[num_vertices*index+i]-1)];
     61                        this->inputs->AddInput(new TriaVertexInput(ControlParameterEnum,nodeinputs));
     62                }
    5763        }
    5864
    … …  
    7480                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->rheology_n[index];
    7581                        this->inputs->AddInput(new PentaVertexInput(RheologyNEnum,nodeinputs));
     82                }
     83
     84                /*Get control_parameter*/
     85                if (iomodel->control_parameter) {
     86                        for(i=0;i<num_vertices;i++)nodeinputs[i]=iomodel->control_parameter[int(iomodel->elements[num_vertices*index+i]-1)];
     87                        this->inputs->AddInput(new PentaVertexInput(ControlParameterEnum,nodeinputs));
    7688                }
    7789        }
    … …  
    531543}
    532544/*}}}*/
     545/*FUNCTION Matice::InputDuplicate{{{1*/
     546void  Matice::InputDuplicate(int original_enum,int new_enum){
     547
     548        /*Call inputs method*/
     549        printf("%s\n",EnumAsString(original_enum));
     550        if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum);
     551
     552}
     553/*}}}*/
    533554/*FUNCTION Matice::InputUpdateFromVector(double* vector, int name, int type) {{{1*/
    534555void  Matice::InputUpdateFromVector(double* vector, int name, int type){
    535         /*Nothing updated yet*/
     556
     557        /*Intermediaries*/
     558        Element *element      = NULL;
     559
     560        /*Recover element*/
     561        element=(Element*)helement->delivers();
     562
     563        /*Check that name is an element input*/
     564        if (!IsInput(name)) return;
     565
     566        switch(type){
     567
     568                case VertexEnum:
     569
     570                        switch(element->Enum()){
     571
     572                                case TriaEnum:
     573                                        double values[3];
     574                                        for (int i=0;i<3;i++) values[i]=vector[((Tria*)element)->nodes[i]->GetVertexDof()];
     575                                        this->inputs->AddInput(new TriaVertexInput(name,values));
     576                                        return;
     577
     578                                default: ISSMERROR("element %s not implemented yet",EnumAsString(element->Enum()));
     579                        }
     580                default: ISSMERROR("type %i (%s) not implemented yet",type,EnumAsString(type));
     581        }
    536582}
    537583/*}}}*/
    … …  
    566612}
    567613/*}}}*/
     614/*FUNCTION Matice::IsInput{{{1*/
     615bool Matice::IsInput(int name){
     616        if (
     617                                name==RheologyBEnum ||
     618                                name==RheologyB2dEnum ||
     619                                name==RheologyNEnum ||
     620                                name==ControlParameterEnum
     621                ){
     622                return true;
     623        }
     624        else return false;
     625}
     626/*}}}*/

  • issm/trunk/src/c/objects/Materials/Matice.h

    r4931 r4967  
    5050                void  InputUpdateFromSolution(double* solution);
    5151                /*}}}*/
     52                /*Material virtual functions resolution: {{{1*/
     53                void   InputDuplicate(int original_enum,int new_enum);
     54                /*}}}*/
    5255                /*Matice Numerics: {{{1*/
    5356                void   Configure(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin);
    … …  
    6063                double GetB2d();
    6164                double GetN();
     65                bool   IsInput(int name);
    6266                /*}}}*/
    6367};

  • issm/trunk/src/c/objects/Materials/Matpar.h

    r4966 r4967  
    5454                void   InputUpdateFromSolution(double* solution);
    5555                /*}}}*/
     56                /*Material virtual functions resolution: {{{1*/
     57                void   InputDuplicate(int original_enum,int new_enum){ISSMERROR("not implemented yet");};
     58                /*}}}*/
    5659                /*Numerics: {{{1*/
    5760                double GetG();

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

    r4894 r4967  
    2323        double  eps_cm;
    2424        double  tolx;
    25         double  cm_min;
    26         double  cm_max;
    2725        bool    cm_gradient;
    2826        int     dim;
    … …  
    5452        femmodel->parameters->FindParam(&eps_cm,EpsCmEnum);
    5553        femmodel->parameters->FindParam(&tolx,TolXEnum);
    56         femmodel->parameters->FindParam(&cm_min,CmMinEnum);
    57         femmodel->parameters->FindParam(&cm_max,CmMaxEnum);
    5854        femmodel->parameters->FindParam(&cm_gradient,CmGradientEnum);
    5955        femmodel->parameters->FindParam(&control_steady,ControlSteadyEnum);
    … …  
    9086                BrentSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs);
    9187
    92                 _printf_("%s\n","      updating parameter using optimized search scalar...");
    93                 InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ControlParameterEnum,control_type);
    94                 InputAXPYx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,search_scalar*optscal[n],GradientEnum);
    95 
    96                 _printf_("%s\n","      constraining the new distribution...");   
    97                 InputControlConstrainx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,cm_min,cm_max);
    98                
    99                 _printf_("%s\n","      save new parameter...");
    100                 InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,ControlParameterEnum);
    101                
    102                 _printf_("%s%i%s%g\n","      value of misfit J after optimization #",n+1,": ",J[n]);
     88                _printf_("%s\n","      updating parameter using optimized search scalar..."); //true means update parameter and copy it onto ControlParameter input
     89                InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar*optscal[n],true);
    10390               
    10491                if(controlconvergence(J,fit,eps_cm,n)) break;

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

    r4839 r4967  
    3232        double* optscal=NULL;
    3333        double* fit=NULL;
    34         double  cm_min;
    35         double  cm_max;
    3634        int     control_type;
    3735        int     analysis_type;
    … …  
    4846        femmodel->parameters->FindParam(&optscal,NULL,OptScalEnum);
    4947        femmodel->parameters->FindParam(&fit,NULL,FitEnum);
    50         femmodel->parameters->FindParam(&cm_min,CmMinEnum);
    51         femmodel->parameters->FindParam(&cm_max,CmMaxEnum);
    5248        femmodel->parameters->FindParam(&isstokes,IsStokesEnum);
    5349        femmodel->parameters->FindParam(&control_type,ControlTypeEnum);
    … …  
    5955        else femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
    6056
    61         /*Use ControlParameterEnum input to  reinitialize our input parameter: */
    62         InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ControlParameterEnum,control_type);
    63        
    64         /*Use search scalar to shoot parameter in the gradient direction: */
    65         InputAXPYx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,search_scalar*optscal[n],GradientEnum);
    66 
    67         /*Constrain:*/
    68         InputControlConstrainx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,cm_min,cm_max);
     57        /*update parameter according to scalar: */ //false means: do not copy updated parameter onto ControlParameter input
     58        InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar*optscal[n],false);
    6959
    7060        /*Run diagnostic with updated inputs: */

  • issm/trunk/src/m/solutions/control_core.m

    r4895 r4967  
    1919        eps_cm=femmodel.parameters.EpsCm;
    2020        tolx=femmodel.parameters.TolX;
    21         cm_min=femmodel.parameters.CmMin;
    22         cm_max=femmodel.parameters.CmMax;
    2321        cm_gradient=femmodel.parameters.CmGradient;
    2422        control_steady=femmodel.parameters.ControlSteady;
    … …  
    5553
    5654                displaystring('\n%s',['      updating parameter using optimized search scalar:']);
    57                 scalar=search_scalar*optscal(n);
    58                 femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ControlParameterEnum,control_type);
    59                 [femmodel.elements femmodel.materials]=InputAXPY(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,scalar,GradientEnum);
    60 
    61                 displaystring('\n%s',['      constraning the new distribution...']);
    62                 [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlConstrain(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,cm_min,cm_max);
    63 
    64                 displaystring('\n%s',['      save new parameter...']);
    65                 femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,ControlParameterEnum);
     55                [femmodel.elements,femmodel.nodes,femmmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlUpdate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,search_scalar*optscal(n),1);
    6656
    6757                disp(['      value of misfit J after optimization #' num2str(n) ':' num2str(J(n))]);

  • issm/trunk/src/m/solutions/objectivefunctionC.m

    r4804 r4967  
    1010analysis_type=femmodel.parameters.AnalysisType;
    1111isstokes=femmodel.parameters.IsStokes;
    12 cm_min=femmodel.parameters.CmMin;
    13 cm_max=femmodel.parameters.CmMax;
    1412
    1513%set current configuration
    … …  
    2018end
    2119
    22 %Use ControlParameterEnum input to  reinitialize our input parameter:
    23 femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ControlParameterEnum,control_type);
    24 
    2520%Use search scalar to shoot parameter in the gradient direction:
    26 scalar=search_scalar*optscal;
    27 [femmodel.elements femmodel.materials]=InputAXPY(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,scalar,GradientEnum);
    28 
    29 %Constrain:
    30 [femmodel.elements,femmodel.nodes,femmmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlConstrain(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,control_type,cm_min,cm_max);
     21[femmodel.elements,femmodel.nodes,femmmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputControlUpdate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,search_scalar*optscal,0);
    3122
    3223%Run diagnostic with updated inputs:

  • issm/trunk/src/mex/Makefile.am

    r4938 r4967  
    2222                                ElementConnectivity\
    2323                                InputControlConstrain \
     24                                InputControlUpdate \
    2425                                InputConvergence\
    2526                                GetSolutionFromInputs\
    … …  
    134135                          InputControlConstrain/InputControlConstrain.h
    135136
     137InputControlUpdate_SOURCES = InputControlUpdate/InputControlUpdate.cpp\
     138                                                                                  InputControlUpdate/InputControlUpdate.h
     139
    136140InputConvergence_SOURCES = InputConvergence/InputConvergence.cpp\
    137141                          InputConvergence/InputConvergence.h
Note: See TracChangeset for help on using the changeset viewer.