Changeset 8646


Ignore:
Timestamp:
06/16/11 14:13:46 (14 years ago)
Author:
Mathieu Morlighem
Message:

Deleted noisedmp, use DragCoefficientAbsGradientEnum as a cost function instead

Location:
issm/trunk/src/c/objects
Files:
6 edited

Legend:

Unmodified
Added
Removed

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

    r8643 r8646  
    38763876/*FUNCTION Penta::Gradj {{{1*/
    38773877void  Penta::Gradj(Vec gradient,int control_type){
     3878        /*dJ/dalpha = ∂L/∂alpha = ∂J/∂alpha + ∂/∂alpha(KU-F)*/
    38783879
    38793880        int              i,approximation;
    … …  
    38833884        if(IsOnWater())return;
    38843885
     3886        /*First deal with ∂/∂alpha(KU-F)*/
    38853887        switch(control_type){
    38863888
    … …  
    39623964                default:
    39633965                        _error_("control type %s not supported yet: ",EnumToStringx(control_type));
     3966        }
     3967
     3968        /*Now deal with ∂J/∂alpha*/
     3969        int        *responses = NULL;
     3970        int         num_responses,resp;
     3971        this->parameters->FindParam(&num_responses,NumResponsesEnum);
     3972        this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
     3973
     3974        for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
     3975
     3976                case ThicknessAbsMisfitEnum:
     3977                case SurfaceAbsVelMisfitEnum:
     3978                case SurfaceRelVelMisfitEnum:
     3979                case SurfaceLogVelMisfitEnum:
     3980                case SurfaceLogVxVyMisfitEnum:
     3981                case SurfaceAverageVelMisfitEnum:
     3982                        /*Nothing, J does not depends on the parameter being inverted for*/
     3983                        break;
     3984                case DragCoefficientAbsGradientEnum:
     3985                        if (!IsOnBed()) return;
     3986                        tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
     3987                        tria->GradjDragGradient(gradient,resp);
     3988                        delete tria->matice; delete tria;
     3989                        break;
     3990                default:
     3991                        _error_("response %s not supported yet",EnumToStringx(responses[resp]));
    39643992        }
    39653993}
    … …  
    50465074                /*Create inputs and add to DataSetInput*/
    50475075                DatasetInput* datasetinput=new DatasetInput(WeightsEnum);
    5048                 for(j=0;j<iomodel->num_cm_responses;j++){
    5049                         for(i=0;i<6;i++)nodeinputs[i]=iomodel->weights[penta_vertex_ids[i]-1];
     5076                for(i=0;i<iomodel->num_cm_responses;i++){
     5077                        for(j=0;j<6;j++)nodeinputs[j]=iomodel->weights[(penta_vertex_ids[j]-1)*iomodel->num_cm_responses+i];
    50505078                        datasetinput->inputs->AddObject(new PentaVertexInput(WeightsEnum,nodeinputs));
    50515079                }

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

    r8643 r8646  
    19901990                                        }
    19911991                                        break;
     1992                                case DragCoefficientAbsGradientEnum:
     1993                                        /*Nothing in P vector*/
     1994                                        break;
     1995                                case ThicknessAbsGradientEnum:
     1996                                        /*Nothing in P vector*/
     1997                                        break;
     1998                                case RheologyBbarAbsGradientEnum:
     1999                                        /*Nothing in P vector*/
     2000                                        break;
    19922001                                default:
    19932002                                        _error_("response %s not supported yet",EnumToStringx(responses[resp]));
    … …  
    28762885/*FUNCTION Tria::Gradj {{{1*/
    28772886void  Tria::Gradj(Vec gradient,int control_type){
     2887        /*dJ/dalpha = ∂L/∂alpha = ∂J/∂alpha + ∂/∂alpha(KU-F)*/
    28782888
    28792889        /*If on water, grad = 0: */
    28802890        if(IsOnWater()) return;
    28812891
     2892        /*First deal with ∂/∂alpha(KU-F)*/
    28822893        switch(control_type){
    28832894                case DragCoefficientEnum:
    … …  
    28992910                        _error_("%s%i","control type not supported yet: ",control_type);
    29002911        }
     2912
     2913        /*Now deal with ∂J/∂alpha*/
     2914        int        *responses = NULL;
     2915        int         num_responses,resp;
     2916        this->parameters->FindParam(&num_responses,NumResponsesEnum);
     2917        this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
     2918
     2919        for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
     2920
     2921                case ThicknessAbsMisfitEnum:
     2922                case SurfaceAbsVelMisfitEnum:
     2923                case SurfaceRelVelMisfitEnum:
     2924                case SurfaceLogVelMisfitEnum:
     2925                case SurfaceLogVxVyMisfitEnum:
     2926                case SurfaceAverageVelMisfitEnum:
     2927                        /*Nothing, J does not depends on the parameter being inverted for*/
     2928                        break;
     2929                case DragCoefficientAbsGradientEnum:
     2930                        GradjDragGradient(gradient,resp);
     2931                        break;
     2932                default:
     2933                        _error_("response %s not supported yet",EnumToStringx(responses[resp]));
     2934        }
    29012935}
    29022936/*}}}*/
    … …  
    31573191
    31583192        delete friction;
     3193        delete gauss;
     3194}
     3195/*}}}*/
     3196/*FUNCTION Tria::GradjDragGradient{{{1*/
     3197void  Tria::GradjDragGradient(Vec gradient, int weight_index){
     3198
     3199        int        i,ig;
     3200        int        doflist1[NUMVERTICES];
     3201        double     Jdet,weight;
     3202        double     xyz_list[NUMVERTICES][3];
     3203        double     dh1dh3[NDOF2][NUMVERTICES];
     3204        double     dk[NDOF2];
     3205        double     grade_g[NUMVERTICES]={0.0};
     3206        GaussTria  *gauss=NULL;
     3207
     3208        /*Retrieve all inputs we will be needing: */
     3209        if(IsOnShelf())return;
     3210        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     3211        GetDofList1(&doflist1[0]);
     3212        Input* dragcoefficient_input=inputs->GetInput(DragCoefficientEnum); _assert_(dragcoefficient_input);
     3213        Input* weights_input=inputs->GetInput(WeightsEnum);                 _assert_(weights_input);
     3214
     3215        /* Start  looping on the number of gaussian points: */
     3216        gauss=new GaussTria(2);
     3217        for (ig=gauss->begin();ig<gauss->end();ig++){
     3218
     3219                gauss->GaussPoint(ig);
     3220
     3221                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     3222                GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss);
     3223                weights_input->GetParameterValue(&weight,gauss,weight_index);
     3224
     3225                /*Build alpha_complement_list: */
     3226                dragcoefficient_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
     3227
     3228                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
     3229                for (i=0;i<NUMVERTICES;i++){
     3230                        grade_g[i]+=-weight*Jdet*gauss->weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
     3231                        _assert_(!isnan(grade_g[i]));
     3232                }
     3233        }
     3234        VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);
     3235
     3236        /*Clean up and return*/
    31593237        delete gauss;
    31603238}
    … …  
    36413719                /*Create inputs and add to DataSetInput*/
    36423720                DatasetInput* datasetinput=new DatasetInput(WeightsEnum);
    3643                 for(j=0;j<iomodel->num_cm_responses;j++){
    3644                         for(i=0;i<3;i++)nodeinputs[i]=iomodel->weights[tria_vertex_ids[i]-1];
     3721                for(i=0;i<iomodel->num_cm_responses;i++){
     3722                        for(j=0;j<3;j++)nodeinputs[j]=iomodel->weights[(tria_vertex_ids[j]-1)*iomodel->num_cm_responses+i];
    36453723                        datasetinput->inputs->AddObject(new TriaVertexInput(WeightsEnum,nodeinputs));
    36463724                }

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

    r8643 r8646  
    9393                void   GradjDragMacAyeal(Vec gradient);
    9494                void   GradjDragStokes(Vec gradient);
     95                void   GradjDragGradient(Vec gradient,int weight_index);
    9596                void   GradjDhDtBalancedthickness(Vec gradient);
    9697                void   GradjVxBalancedthickness(Vec gradient);

  • issm/trunk/src/c/objects/Params/DoubleMatParam.cpp

    r8224 r8646  
    152152
    153153/*DoubleMatParam virtual functions definitions: */
    154 /*FUNCTION DoubleMatParam::GetParameterValue{{{1*/
     154/*FUNCTION DoubleMatParam::GetParameterValue(double** pdoublearray,int* pM,int* pN){{{1*/
    155155void  DoubleMatParam::GetParameterValue(double** pdoublearray,int* pM,int* pN){
    156156        double* output=NULL;
    … …  
    163163        if(pN) *pN=N;
    164164        *pdoublearray=output;
     165}
     166/*}}}*/
     167/*FUNCTION DoubleMatParam::GetParameterValue(int** pintarray,int* pM,int* pN){{{1*/
     168void  DoubleMatParam::GetParameterValue(int** pintarray,int* pM,int* pN){
     169#ifdef _SERIAL_
     170        int* output=NULL;
     171        int  i;
     172
     173        output=(int*)xmalloc((int)(M*N*sizeof(int)));
     174        for(i=0;i<M*N;i++) output[i]=(int)value[i];
     175
     176        /*Assign output pointers:*/
     177        if(pM) *pM=M;
     178        if(pN) *pN=N;
     179        *pintarray=output;
     180#else
     181        _error_("DoubleMat of enum %i (%s) cannot return an array of int",enum_type,EnumToStringx(enum_type));
     182#endif
    165183}
    166184/*}}}*/

  • issm/trunk/src/c/objects/Params/DoubleMatParam.h

    r8600 r8646  
    5454                void  GetParameterValue(int* pinteger){_error_("DoubleMat param of enum %i (%s) cannot return an integer",enum_type,EnumToStringx(enum_type));}
    5555                void  GetParameterValue(int** pintarray,int* pM){_error_("DoubleMat param of enum %i (%s) cannot return an array of integers",enum_type,EnumToStringx(enum_type));}
    56                 void  GetParameterValue(int** pintarray,int* pM,int* pN){_error_("DoubleMat param of enum %i (%s) cannot return an array of integers",enum_type,EnumToStringx(enum_type));}
     56                void  GetParameterValue(int** pintarray,int* pM,int* pN);
    5757                void  GetParameterValue(double* pdouble){_error_("DoubleMat param of enum %i (%s) cannot return a double",enum_type,EnumToStringx(enum_type));}
    5858                void  GetParameterValue(char** pstring){_error_("DoubleMat param of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));}

  • issm/trunk/src/c/objects/Params/DoubleVecParam.cpp

    r8224 r8646  
    174174        *pintarray=output;
    175175#else
    176         _error_("Double param of enum %i (%s) cannot return an array of double",enum_type,EnumToStringx(enum_type));
     176        _error_("DoubleVec param of enum %i (%s) cannot return an array of int",enum_type,EnumToStringx(enum_type));
    177177#endif
    178178}
Note: See TracChangeset for help on using the changeset viewer.