Changeset 11317


Ignore:
Timestamp:
02/03/12 13:50:06 (13 years ago)
Author:
Mathieu Morlighem
Message:

Added 2 modules to make TAO work

Location:
issm/trunk-jpl/src/c
Files:
6 added
13 edited

Legend:

Unmodified
Added
Removed

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

    r11295 r11317  
    423423                                          ./modules/ControlInputScaleGradientx/ControlInputScaleGradientx.cpp\
    424424                                          ./modules/ControlInputScaleGradientx/ControlInputScaleGradientx.h\
     425                                          ./modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.cpp\
     426                                          ./modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.h\
     427                                          ./modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp\
     428                                          ./modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.h\
    425429                                          ./modules/ModelProcessorx/Control/CreateParametersControl.cpp\
    426430                                          ./modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp\

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

    r8387 r11317  
    1010
    1111void GetVectorFromInputsx( Vec* pvector, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters, int name, int type){
    12 
    1312
    1413        int i;
    … …  
    4140        /*Assign output pointers:*/
    4241        *pvector=vector;
    43 
    4442}
    4543

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

    r11313 r11317  
    2727        /*Allocate gradient_list */
    2828        gradient_list = (Vec*)xmalloc(num_controls*sizeof(Vec));
    29         norm_list     = (double*)xmalloc(num_controls*sizeof(double));
     29        norm_list = (double*)xmalloc(num_controls*sizeof(double));
    3030        for(i=0;i<num_controls;i++){
    3131                gradient_list[i]=NewVec(num_controls*numberofvertices);
    … …  
    5959
    6060        /*Clean-up and assign output pointer*/
    61         *pnorm_list=norm_list;
    62         *pgradient=gradient;
     61        if(pnorm_list){
     62                *pnorm_list=norm_list;
     63        }
     64        else{
     65                xfree((void**)&norm_list);
     66        }
     67        if(pgradient)  *pgradient=gradient;
    6368        xfree((void**)&gradient_list);
    6469        xfree((void**)&control_type);

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

    r11295 r11317  
    3131#include "./GetSolutionFromInputsx/GetSolutionFromInputsx.h"
    3232#include "./GetVectorFromInputsx/GetVectorFromInputsx.h"
     33#include "./GetVectorFromControlInputsx/GetVectorFromControlInputsx.h"
     34#include "./SetControlInputsFromVectorx/SetControlInputsFromVectorx.h"
    3335#include "./Gradjx/Gradjx.h"
    3436#include "./GroundinglineMigrationx/GroundinglineMigrationx.h"

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

    r11313 r11317  
    103103                virtual void   ControlInputSetGradient(double* gradient,int enum_type,int control_index)=0;
    104104                virtual void   ControlInputScaleGradient(int enum_type, double scale)=0;
     105                virtual void   GetVectorFromControlInputs(Vec gradient,int control_enum,int control_index)=0;
     106                virtual void   SetControlInputsFromVector(double* vector,int control_enum,int control_index)=0;
    105107                virtual void   InputControlUpdate(double scalar,bool save_parameter)=0;
    106108                #endif

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

    r11314 r11317  
    16071607        this->results->AddObject((Object*)input->SpawnResult(step,time));
    16081608        #ifdef _HAVE_CONTROL_
    1609         if(input->ObjectEnum()==ControlInputEnum) this->results->AddObject((Object*)((ControlInput*)input)->SpawnGradient(step,time));
     1609        if(input->ObjectEnum()==ControlInputEnum){
     1610                if(((ControlInput*)input)->gradient!=NULL) this->results->AddObject((Object*)((ControlInput*)input)->SpawnGradient(step,time));
     1611        }
    16101612        #endif
    16111613
    … …  
    51295131}
    51305132/*}}}*/
     5133/*FUNCTION Penta::GetVectorFromControlInputs{{{1*/
     5134void  Penta::GetVectorFromControlInputs(Vec vector,int control_enum,int control_index){
     5135
     5136        int doflist1[NUMVERTICES];
     5137
     5138        /*Get out if this is not an element input*/
     5139        if(!IsInput(control_enum)) return;
     5140
     5141        /*Prepare index list*/
     5142        GradientIndexing(&doflist1[0],control_index);
     5143
     5144        /*Get input (either in element or material)*/
     5145        Input* input=inputs->GetInput(control_enum);
     5146        if(!input) _error_("Input %s not found in element",EnumToStringx(control_enum));
     5147
     5148        /*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
     5149        input->GetVectorFromInputs(vector,&doflist1[0]);
     5150}
     5151/*}}}*/
     5152/*FUNCTION Penta::SetControlInputsFromVector{{{1*/
     5153void  Penta::SetControlInputsFromVector(double* vector,int control_enum,int control_index){
     5154
     5155        double  values[NUMVERTICES];
     5156        int     doflist1[NUMVERTICES];
     5157        Input  *input     = NULL;
     5158        Input  *new_input = NULL;
     5159
     5160        /*Get out if this is not an element input*/
     5161        if(!IsInput(control_enum)) return;
     5162
     5163        /*Prepare index list*/
     5164        GradientIndexing(&doflist1[0],control_index);
     5165
     5166        /*Get values on vertices*/
     5167        for (int i=0;i<NUMVERTICES;i++){
     5168                values[i]=vector[doflist1[i]];
     5169        }
     5170        new_input = new PentaP1Input(control_enum,values);
     5171
     5172
     5173        if(control_enum==MaterialsRheologyBbarEnum){
     5174                input=(Input*)matice->inputs->GetInput(control_enum); _assert_(input);
     5175        }
     5176        else{
     5177                input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
     5178        }
     5179
     5180        if (input->ObjectEnum()!=ControlInputEnum){
     5181                _error_("input %s is not a ControlInput",EnumToStringx(control_enum));
     5182        }
     5183
     5184        ((ControlInput*)input)->SetInput(new_input);
     5185}
     5186/*}}}*/
    51315187#endif
    51325188

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

    r11313 r11317  
    149149                void   GradjBbarPattyn(Vec gradient,int control_index);
    150150                void   GradjBbarStokes(Vec gradient,int control_index);
     151                void   GetVectorFromControlInputs(Vec gradient,int control_enum,int control_index);
     152                void   SetControlInputsFromVector(double* vector,int control_enum,int control_index);
    151153                void   ControlInputGetGradient(Vec gradient,int enum_type,int control_index);
    152154                void   ControlInputScaleGradient(int enum_type,double scale);

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

    r11313 r11317  
    14501450       
    14511451        #ifdef _HAVE_CONTROL_
    1452         if(input->ObjectEnum()==ControlInputEnum) this->results->AddObject((Object*)((ControlInput*)input)->SpawnGradient(step,time));
     1452        if(input->ObjectEnum()==ControlInputEnum){
     1453                if(((ControlInput*)input)->gradient!=NULL) this->results->AddObject((Object*)((ControlInput*)input)->SpawnGradient(step,time));
     1454        }
    14531455        #endif
    14541456}
    … …  
    47284730}
    47294731/*}}}*/
     4732/*FUNCTION Tria::GetVectorFromControlInputs{{{1*/
     4733void  Tria::GetVectorFromControlInputs(Vec vector,int control_enum,int control_index){
     4734
     4735        int doflist1[NUMVERTICES];
     4736
     4737        /*Get out if this is not an element input*/
     4738        if(!IsInput(control_enum)) return;
     4739
     4740        /*Prepare index list*/
     4741        GradientIndexing(&doflist1[0],control_index);
     4742
     4743        /*Get input (either in element or material)*/
     4744        Input* input=inputs->GetInput(control_enum);
     4745        if(!input) _error_("Input %s not found in element",EnumToStringx(control_enum));
     4746
     4747        /*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
     4748        input->GetVectorFromInputs(vector,&doflist1[0]);
     4749}
     4750/*}}}*/
     4751/*FUNCTION Tria::SetControlInputsFromVector{{{1*/
     4752void  Tria::SetControlInputsFromVector(double* vector,int control_enum,int control_index){
     4753
     4754        double  values[NUMVERTICES];
     4755        int     doflist1[NUMVERTICES];
     4756        Input  *input     = NULL;
     4757        Input  *new_input = NULL;
     4758
     4759        /*Get out if this is not an element input*/
     4760        if(!IsInput(control_enum)) return;
     4761
     4762        /*Prepare index list*/
     4763        GradientIndexing(&doflist1[0],control_index);
     4764
     4765        /*Get values on vertices*/
     4766        for (int i=0;i<NUMVERTICES;i++){
     4767                values[i]=vector[doflist1[i]];
     4768        }
     4769        new_input = new TriaP1Input(control_enum,values);
     4770
     4771
     4772        if(control_enum==MaterialsRheologyBbarEnum){
     4773                input=(Input*)matice->inputs->GetInput(control_enum); _assert_(input);
     4774        }
     4775        else{
     4776                input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
     4777        }
     4778
     4779        if (input->ObjectEnum()!=ControlInputEnum){
     4780                _error_("input %s is not a ControlInput",EnumToStringx(control_enum));
     4781        }
     4782
     4783        ((ControlInput*)input)->SetInput(new_input);
     4784}
     4785/*}}}*/
    47304786#endif
    47314787

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

    r11313 r11317  
    150150                void   GradjVxBalancedthickness(Vec gradient,int control_index);
    151151                void   GradjVyBalancedthickness(Vec gradient,int control_index);
     152                void   GetVectorFromControlInputs(Vec gradient,int control_enum,int control_index);
     153                void   SetControlInputsFromVector(double* vector,int control_enum,int control_index);
    152154                void   ControlInputGetGradient(Vec gradient,int enum_type,int control_index);
    153155                void   ControlInputScaleGradient(int enum_type,double scale);

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

    r11291 r11317  
    386386
    387387}/*}}}*/
     388/*FUNCTION ControlInput::SetInput{{{1*/
     389void ControlInput::SetInput(Input* in_input){
     390
     391        delete values; this->values=in_input;
     392        this->SaveValue(); //because this is what SpawnResult saves FIXME
     393
     394}/*}}}*/
    388395/*FUNCTION ControlInput::SpawnResult{{{1*/
    389396ElementResult* ControlInput::SpawnResult(int step, double time){
    … …  
    396403/*FUNCTION ControlInput::SpawnGradient{{{1*/
    397404ElementResult* ControlInput::SpawnGradient(int step, double time){
     405        _assert_(gradient);
    398406        return gradient->SpawnResult(step,time);
    399407}/*}}}*/

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

    r10135 r11317  
    5454                /*}}}*/
    5555                /*numerics: {{{1*/
     56                void SetInput(Input* in_input);
    5657                void GetInputValue(bool* pvalue);
    5758                void GetInputValue(int* pvalue);

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

    r11306 r11317  
    7070        /*Initialize responses: */
    7171        J=(double*)xmalloc(nsteps*sizeof(double));
    72         step_responses=(int*)xmalloc(num_responses*sizeof(double));
     72        step_responses=(int*)xmalloc(num_responses*sizeof(int));
    7373               
    7474        /*Initialize some of the BrentSearch arguments: */

  • issm/trunk-jpl/src/c/solutions/controltao_core.cpp

    r11313 r11317  
    2525
    2626        /*TAO*/
    27         int       ierr,numberofvertices;
    28         AppCtx    user;
    29         TaoSolver tao;
    30         Vec       X  = NULL;
    31         Vec       XL = NULL;
    32         Vec       XU = NULL;
     27        int        ierr,numberofvertices;
     28        int        num_controls;
     29        AppCtx     user;
     30        TaoSolver  tao;
     31        int       *control_list   = NULL;
     32        Vec        X              = NULL;
     33        Vec        XL             = NULL;
     34        Vec        XU             = NULL;
    3335
    3436        /*Initialize TAO*/
    … …  
    3739        ierr = TaoInitialize(&argc,&args,(char*)0,"");
    3840        if(ierr) _error_("Could not initialize Tao");
     41
     42        /*Recover some parameters*/
     43        femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     44        femmodel->parameters->FindParam(&control_list,NULL,InversionControlParametersEnum);
    3945
    4046        /*Initialize TAO*/
    … …  
    5258        XL=NewVec(numberofvertices); VecSet(XL,1.);
    5359        XU=NewVec(numberofvertices); VecSet(XU,200.);
    54         TaoSetVariableBounds(tao,XL,XU); VecFree(&XL); VecFree(&XU);
     60        //TaoSetVariableBounds(tao,XL,XU); VecFree(&XL); VecFree(&XU);
    5561
    56         GetVectorFromInputsx(&X,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,FrictionCoefficientEnum,VertexEnum);
     62        GetVectorFromControlInputsx(&X,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    5763        TaoSetInitialVector(tao,X);
    5864
    … …  
    6470        TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);
    6571        TaoGetSolutionVector(tao,&X);
    66         InputUpdateFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X,FrictionCoefficientEnum,VertexEnum);
    67         InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,FrictionCoefficientEnum);
     72        SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
     73        for(int i=0;i<num_controls;i++){
     74                InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_list[i]);
     75        }
    6876
    6977        /*Clean up and return*/
     78        xfree((void**)&control_list);
    7079        VecFree(&X);
    7180        TaoDestroy(&tao);
    … …  
    7584
    7685        /*Retreive arguments*/
    77         AppCtx   *user     = (AppCtx*)userCtx;
    78         FemModel *femmodel = user->femmodel;
    79         Vec       gradient = NULL;
     86        int       solution_type,num_cost_functions;
     87        AppCtx   *user           = (AppCtx *)userCtx;
     88        FemModel *femmodel       = user->femmodel;
     89        int      *cost_functions = NULL;
     90        double   *cost_functionsd= NULL;
     91        Vec       gradient       = NULL;
    8092
    81         int costfunction=SurfaceAbsVelMisfitEnum;
    82         femmodel->parameters->SetParam(&costfunction,1,1,StepResponsesEnum);
    83         //InputUpdateFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X,MaterialsRheologyBbarEnum,VertexEnum);
    84         InputUpdateFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X,FrictionCoefficientEnum,VertexEnum);
    85         adjointdiagnostic_core(user->femmodel);
    86         //Gradjx(&gradient,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MaterialsRheologyBbarEnum);
     93        /*Set new variable*/
     94        SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
     95
     96        /*Recover some parameters*/
     97        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
     98        femmodel->parameters->FindParam(&num_cost_functions,InversionNumCostFunctionsEnum);
     99        femmodel->parameters->FindParam(&cost_functionsd,NULL,NULL,InversionCostFunctionsEnum);
     100
     101        /*Prepare objective function*/
     102        cost_functions=(int*)xmalloc(num_cost_functions*sizeof(int));
     103        for(int i=0;i<num_cost_functions;i++) cost_functions[i]=(int)cost_functionsd[i]; //FIXME
     104        femmodel->parameters->SetParam(cost_functions,1,num_cost_functions,StepResponsesEnum);
     105
     106        /*Compute solution and adjoint*/
     107        void (*solutioncore)(FemModel*)=NULL;
     108        void (*adjointcore)(FemModel*)=NULL;
     109        CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
     110        AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type);
     111        solutioncore(femmodel);
     112        adjointcore(femmodel);
     113
     114        /*Compute objective function*/
     115        CostFunctionx(fcn,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     116
     117        /*Compute gradient*/
    87118        Gradjx(&gradient,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    88         VecScale(gradient,-1.);
    89         VecCopy(gradient,G);
    90         VecFree(&gradient);
    91         CostFunctionx(fcn,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     119        VecCopy(gradient,G); VecFree(&gradient);
     120        VecScale(G,-1.);
     121
     122        /*Clean-up and return*/
     123        xfree((void**)&cost_functions);
     124        xfree((void**)&cost_functionsd);
    92125        return 0;
    93126}
Note: See TracChangeset for help on using the changeset viewer.