Changeset 18565


Ignore:
Timestamp:
10/03/14 08:36:31 (10 years ago)
Author:
Mathieu Morlighem
Message:

CHG: working on new implementation of balancethickness2

Location:
issm/trunk-jpl/src
Files:
4 added
15 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp

    r18476 r18565  
    4545}/*}}}*/
    4646ElementVector* AdjointBalancethickness2Analysis::CreatePVector(Element* element){/*{{{*/
    47         _error_("Not implemented yet");
     47
     48        /*Intermediaries */
     49        int         num_responses,i;
     50        IssmDouble  vx,vy,vel,Jdet;
     51        IssmDouble  surface,surfaceobs,weight;
     52        int        *responses = NULL;
     53        IssmDouble *xyz_list  = NULL;
     54
     55        /*Fetch number of nodes and dof for this finite element*/
     56        int numnodes = element->GetNumberOfNodes();
     57
     58        /*Initialize Element vector and vectors*/
     59        ElementVector* pe     = element->NewElementVector(SSAApproximationEnum);
     60        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     61        IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
     62
     63        /*Retrieve all inputs and parameters*/
     64        element->GetVerticesCoordinates(&xyz_list);
     65        element->FindParam(&num_responses,InversionNumCostFunctionsEnum);
     66        element->FindParam(&responses,NULL,InversionCostFunctionsEnum);
     67        Input* surface_input      = element->GetInput(SurfaceEnum);                          _assert_(surface_input);
     68        Input* surfaceobs_input   = element->GetInput(InversionSurfaceObsEnum);              _assert_(surfaceobs_input);
     69        Input* weights_input      = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     70        Input* vx_input           = element->GetInput(VxEnum);                                 _assert_(vx_input);
     71        Input* vy_input           = element->GetInput(VyEnum);                                 _assert_(vy_input);
     72
     73        /* Start  looping on the number of gaussian points: */
     74        Gauss* gauss=element->NewGauss(2);
     75        for(int ig=gauss->begin();ig<gauss->end();ig++){
     76                gauss->GaussPoint(ig);
     77
     78                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     79                element->NodalFunctions(basis,gauss);
     80                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     81
     82                surface_input->GetInputValue(&surface, gauss);
     83                surfaceobs_input->GetInputValue(&surfaceobs, gauss);
     84
     85                /*Loop over all requested responses*/
     86                for(int resp=0;resp<num_responses;resp++){
     87                        weights_input->GetInputValue(&weight,gauss,responses[resp]);
     88
     89                        switch(responses[resp]){
     90                                case SurfaceAbsMisfitEnum:
     91                                        for(i=0;i<numnodes;i++) pe->values[i]+=(surfaceobs-surface)*weight*Jdet*gauss->weight*basis[i];
     92                                        break;
     93                                default:
     94                                        _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
     95                        }
     96                }
     97        }
     98
     99        /*Clean up and return*/
     100        xDelete<int>(responses);
     101        xDelete<IssmDouble>(xyz_list);
     102        xDelete<IssmDouble>(basis);
     103        xDelete<IssmDouble>(dbasis);
     104        delete gauss;
     105        return pe;
    48106
    49107}/*}}}*/
    … …  
    71129        /*Deal with first part (partial derivative a J with respect to k)*/
    72130        for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
    73                 case Balancethickness2MisfitEnum:
     131                case SurfaceAbsMisfitEnum:
    74132                        /*Nothing, \partial J/\partial k = 0*/
    75133                        break;
    … …  
    79137        /*Deal with second term*/
    80138        switch(control_type){
     139                case BalancethicknessOmegaEnum:  GradientJOmega(element,gradient,control_index); break;
    81140                default: _error_("control type not supported yet: " << EnumToStringx(control_type));
    82141        }
    … …  
    86145
    87146}/*}}}*/
    88 void AdjointBalancethickness2Analysis::GradientJAdot(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
     147void AdjointBalancethickness2Analysis::GradientJOmega(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    89148
    90149        /*Intermediaries*/
    91         IssmDouble lambda,Jdet;
     150        IssmDouble dlambda[2],ds[2],D0,Jdet;
    92151        IssmDouble *xyz_list= NULL;
    93152
    … …  
    103162        element->GetVerticesCoordinates(&xyz_list);
    104163        element->GradientIndexing(&vertexpidlist[0],control_index);
    105         Input* adjoint_input = element->GetInput(AdjointEnum); _assert_(adjoint_input);
     164        Input* adjoint_input = element->GetInput(AdjointEnum);            _assert_(adjoint_input);
     165        Input* s_input       = element->GetInput(SurfaceEnum);            _assert_(s_input);
     166        Input* D0_input      = element->GetInput(BalancethicknessD0Enum); _assert_(D0_input);
    106167
    107168        Gauss* gauss=element->NewGauss(2);
    … …  
    111172                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    112173                element->NodalFunctionsP1(basis,gauss);
    113                 adjoint_input->GetInputValue(&lambda,gauss);
     174
     175                D0_input->GetInputValue(&D0,gauss);
     176                adjoint_input->GetInputDerivativeValue(&dlambda[0],xyz_list,gauss);
     177                s_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss);
    114178
    115179                /*Build gradient vector (actually -dJ/da): */
    116180                for(int i=0;i<numvertices;i++){
    117                         ge[i]+= - Jdet*gauss->weight*basis[i]*lambda;
     181                        ge[i]+= Jdet*gauss->weight*basis[i]*(ds[0]*dlambda[0] + ds[1]*dlambda[1]);
    118182                        _assert_(!xIsNan<IssmDouble>(ge[i]));
    119183                }

  • issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.h

    r18476 r18565  
    2828                void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    2929                void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
    30                 void GradientJAdot(Element* element,Vector<IssmDouble>* gradient,int control_index);
     30                void GradientJOmega(Element* element,Vector<IssmDouble>* gradient,int control_index);
    3131                void InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3232                void UpdateConstraints(FemModel* femmodel);

  • issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp

    r18551 r18565  
    1515        /*Finite element type*/
    1616        int finiteelement = P1Enum;
     17
     18        /*Load variables in element*/
     19        iomodel->FetchDataToInput(elements,ThicknessEnum);
     20        iomodel->FetchDataToInput(elements,SurfaceEnum);
     21        iomodel->FetchDataToInput(elements,BaseEnum);
     22        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
     23        iomodel->FetchDataToInput(elements,BasalforcingsGroundediceMeltingRateEnum);
     24        iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum);
     25        iomodel->FetchDataToInput(elements,BalancethicknessThickeningRateEnum);
     26        iomodel->FetchDataToInput(elements,BalancethicknessOmegaEnum);
    1727
    1828        /*Update elements: */
    … …  
    2232                        Element* element=(Element*)elements->GetObjectByOffset(counter);
    2333                        element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement);
     34
    2435                        counter++;
    2536                }
    2637        }
    2738
    28         iomodel->FetchDataToInput(elements,ThicknessEnum);
    29         iomodel->FetchDataToInput(elements,SurfaceEnum);
    30         iomodel->FetchDataToInput(elements,BaseEnum);
    31         iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
    32         iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
    33         iomodel->FetchDataToInput(elements,BasalforcingsGroundediceMeltingRateEnum);
    34         iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum);
    35         iomodel->FetchDataToInput(elements,BalancethicknessThickeningRateEnum);
    36         iomodel->FetchDataToInput(elements,BalancethicknessCmuEnum);
    3739}/*}}}*/
    3840void Balancethickness2Analysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
    … …  
    6668
    6769        /*Intermediaries */
    68         IssmDouble  Jdet,D;
     70        IssmDouble  Jdet,D0,omega;
    6971        IssmDouble* xyz_list = NULL;
    7072
    … …  
    7678        IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
    7779
    78         /*Create input D*/
    79         this->CreateDiffusionCoefficient(element);
    80 
    8180        /*Retrieve all inputs and parameters*/
    8281        element->GetVerticesCoordinates(&xyz_list);
    83         Input* D_input = element->GetInput(BalancethicknessDiffusionCoefficientEnum); _assert_(D_input);
     82        Input* omega_input = element->GetInput(BalancethicknessOmegaEnum); _assert_(omega_input);
     83        Input* D0_input    = element->GetInput(BalancethicknessD0Enum);
     84        if(!D0_input){
     85                this->CreateD0(element);
     86                D0_input = element->GetInput(BalancethicknessD0Enum); _assert_(D0_input);
     87        }
    8488
    8589        /* Start  looping on the number of gaussian points: */
    … …  
    8993                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    9094                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    91                 D_input->GetInputValue(&D,gauss);
     95                D0_input->GetInputValue(&D0,gauss);
     96                omega_input->GetInputValue(&omega,gauss);
    9297
    9398                for(int i=0;i<numnodes;i++){
    9499                        for(int j=0;j<numnodes;j++){
    95                                 Ke->values[i*numnodes+j] += D*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]);
     100                                Ke->values[i*numnodes+j] += D0*omega*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]);
    96101                        }
    97102                }
    … …  
    107112
    108113        /*Intermediaries */
    109         IssmDouble  dhdt,mb,ms,D,db[2],Jdet;
     114        IssmDouble  dhdt,mb,ms,Jdet;
    110115        IssmDouble* xyz_list = NULL;
    111116
    … …  
    116121        ElementVector* pe    = element->NewElementVector();
    117122        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    118         IssmDouble*   dbasis = xNew<IssmDouble>(2*numnodes);
    119123
    120124        /*Retrieve all inputs and parameters*/
    … …  
    123127        Input* mb_input   = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);       _assert_(mb_input);
    124128        Input* dhdt_input = element->GetInput(BalancethicknessThickeningRateEnum);            _assert_(dhdt_input);
    125         Input* D_input   = element->GetInput(BalancethicknessDiffusionCoefficientEnum);       _assert_(D_input);
    126         Input* bed_input  = element->GetInput(BaseEnum);                                      _assert_(bed_input);
    127129
    128130        /* Start  looping on the number of gaussian points: */
    … …  
    133135                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    134136                element->NodalFunctions(basis,gauss);
    135                 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    136137
    137138                ms_input->GetInputValue(&ms,gauss);
    138139                mb_input->GetInputValue(&mb,gauss);
    139140                dhdt_input->GetInputValue(&dhdt,gauss);
    140                 D_input->GetInputValue(&D,gauss);
    141                 bed_input->GetInputDerivativeValue(&db[0],xyz_list,gauss);
    142                 //db[0]=0.;
    143                 //db[1]=0.;
    144141
    145142                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(
    146143                                        (ms-mb-dhdt)*basis[i]
    147                                         //-D*db[0]*dbasis[0*numnodes+i] - D*db[1]*dbasis[1*numnodes+i]
    148144                                        );
    149145        }
    … …  
    220216
    221217/*Specifics*/
    222 void Balancethickness2Analysis::CreateDiffusionCoefficient(Element* element){/*{{{*/
     218void Balancethickness2Analysis::CreateD0(Element* element){/*{{{*/
    223219
    224220        /*Intermediaries */
    … …  
    231227        /*Fetch number of vertices and allocate output*/
    232228        int  numnodes = element->GetNumberOfNodes();
    233         IssmDouble* D      = xNew<IssmDouble>(numnodes);
    234         IssmDouble* Dgradb = xNew<IssmDouble>(numnodes);
     229        IssmDouble* D0     = xNew<IssmDouble>(numnodes);
    235230
    236231        /*retrieve what we need: */
    237232        element->GetVerticesCoordinates(&xyz_list);
    238         Input* bed_input        = element->GetInput(BaseEnum);               _assert_(bed_input);
    239233        Input* surface_input    = element->GetInput(SurfaceEnum);            _assert_(surface_input);
    240234        Input* B_input          = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
    241         Input* k_input          = element->GetInput(FrictionCoefficientEnum);_assert_(k_input);
    242         Input* Cmu_input        = element->GetInput(BalancethicknessCmuEnum);_assert_(Cmu_input);
    243235        IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
    244236
    … …  
    249241               
    250242                B_input->GetInputValue(&B,gauss);
    251                 k_input->GetInputValue(&k,gauss);
    252                 //thickness_input->GetInputValue(&h,gauss);
    253243                surface_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss);
    254                 surface_input->GetInputValue(&s,gauss);
    255                 bed_input->GetInputValue(&b,gauss);
    256                 h = s-b;
    257244
    258245                mu0   = pow(2.,(1-3*n)/(2.*n)) * B;
    259246                Gamma = pow(rhog,n) * pow(Hstar,2*(n+1)) * pow(Hstar/Lstar,2*(n+1)) * 1./pow(mu0,n);
    260                 //Cmu   = pow(mu0,n)/Hstar/(k*k);
    261                 Cmu_input->GetInputValue(&Cmu,gauss);
    262 
    263                 D[i] = Gamma*pow(ds[0]*ds[0]+ds[1]*ds[1],(n-1)/2)*pow(h,n+1)*(1./(n+2.)*h + Cmu);
     247
     248                D0[i] = Gamma*pow(ds[0]*ds[0]+ds[1]*ds[1],(n-1)/2)/(n+2);
    264249        }
    265250
    266251        /*Add input*/
    267         element->AddInput(BalancethicknessDiffusionCoefficientEnum,D,element->GetElementType());
     252        element->AddInput(BalancethicknessD0Enum,D0,element->GetElementType());
    268253       
    269254        /*Clean up and return*/
    270         xDelete<IssmDouble>(D);
     255        xDelete<IssmDouble>(D0);
    271256        delete gauss;
    272257}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.h

    r18476 r18565  
    3232
    3333                /*Specifics*/
    34                 void CreateDiffusionCoefficient(Element* element);
     34                void CreateD0(Element* element);
    3535};
    3636#endif

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

    r18534 r18565  
    10341034                                /*No yts conversion*/
    10351035                                case ThicknessEnum:
     1036                                case BalancethicknessOmegaEnum:
    10361037                                case FrictionCoefficientEnum:
    10371038                                        if(iomodel->Data(control)){

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

    r18530 r18565  
    641641                case DragCoefficientAbsGradientEnum:DragCoefficientAbsGradientx(responses, elements,nodes, vertices, loads, materials, parameters); break;
    642642                case BalancethicknessMisfitEnum:   BalancethicknessMisfitx(responses); break;
    643                 case Balancethickness2MisfitEnum:  Balancethickness2Misfitx(responses); break;
    644643                case TotalSmbEnum:                                        this->TotalSmbx(responses); break;
    645644                case MaterialsRheologyBbarEnum:    this->ElementResponsex(responses,MaterialsRheologyBbarEnum); break;
    … …  
    719718                                case DragCoefficientAbsGradientEnum:DragCoefficientAbsGradientx(&double_result,elements,nodes,vertices,loads,materials,parameters); break;
    720719                                case BalancethicknessMisfitEnum:    BalancethicknessMisfitx(&double_result);                                                        break;
    721                                 case Balancethickness2MisfitEnum:   Balancethickness2Misfitx(&double_result); break;
     720                                case SurfaceAbsMisfitEnum:          SurfaceAbsMisfitx(&double_result); break;
    722721
    723722                           /*Vector */
    … …  
    14441443
    14451444}/*}}}*/
    1446 void FemModel::Balancethickness2Misfitx(IssmDouble* presponse){/*{{{*/
     1445void FemModel::SurfaceAbsMisfitx(IssmDouble* presponse){/*{{{*/
    14471446
    14481447        /*output: */
    … …  
    14501449        IssmDouble J_sum;
    14511450
    1452         IssmDouble  weight,thicknessobs,thickness,potential,dpotential[2];
    1453         IssmDouble  vx,vy,vxbar,vybar,vxobs,vyobs,vxbarobs,vybarobs,nu;
     1451        IssmDouble  surface,surfaceobs,weight;
    14541452        IssmDouble  Jdet;
    14551453        IssmDouble* xyz_list = NULL;
    … …  
    14611459                /*If on water, return 0: */
    14621460                if(!element->IsIceInElement()) continue;
    1463                  _error_("Not implemented");
     1461
     1462                /* Get node coordinates*/
     1463                element->GetVerticesCoordinates(&xyz_list);
     1464
     1465                 /*Retrieve all inputs we will be needing: */
     1466                 Input* weights_input   =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     1467                 Input* surface_input   =element->GetInput(SurfaceEnum);                            _assert_(surface_input);
     1468                 Input* surfaceobs_input=element->GetInput(InversionSurfaceObsEnum);                _assert_(surfaceobs_input);
     1469
     1470                 /* Start  looping on the number of gaussian points: */
     1471                 Gauss* gauss=element->NewGauss(2);
     1472                 for(int ig=gauss->begin();ig<gauss->end();ig++){
     1473
     1474                         gauss->GaussPoint(ig);
     1475
     1476                         /* Get Jacobian determinant: */
     1477                         element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     1478
     1479                         /*Get all parameters at gaussian point*/
     1480                         weights_input->GetInputValue(&weight,gauss,SurfaceAbsMisfitEnum);
     1481                         surface_input->GetInputValue(&surface,gauss);
     1482                         surfaceobs_input->GetInputValue(&surfaceobs,gauss);
     1483
     1484                         /*Compute SurfaceAbsMisfitEnum*/
     1485                         J+=0.5*(surface-surfaceobs)*(surface-surfaceobs)*weight*Jdet*gauss->weight;
     1486                 }
    14641487
    14651488        }

  • issm/trunk-jpl/src/c/classes/FemModel.h

    r18529 r18565  
    8181                void ElementResponsex(IssmDouble* presponse,int response_enum);
    8282                void BalancethicknessMisfitx(IssmDouble* pV);
    83                 void Balancethickness2Misfitx(IssmDouble* pV);
    8483                #ifdef  _HAVE_DAKOTA_
    8584                void DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses);
    … …  
    9392                void CostFunctionx(IssmDouble* pJ,IssmDouble** pJlist,int* pn);
    9493                void ThicknessAbsGradientx( IssmDouble* pJ);
     94                void SurfaceAbsMisfitx( IssmDouble* pJ);
    9595                #ifdef _HAVE_GIA_
    9696                void Deflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, IssmDouble* x, IssmDouble* y);

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

    r18551 r18565  
    2222
    2323        if(VerboseSolution()) _printf0_("call computational core:\n");
    24         //solutionsequence_linear(femmodel);
    25         solutionsequence_nonlinear(femmodel,false);
     24        solutionsequence_linear(femmodel);
     25        //solutionsequence_nonlinear(femmodel,false);
    2626
    2727        if(save_results){
    2828                if(VerboseSolution()) _printf0_("   saving results\n");
    29                 const int numoutputs = 5;
    30                 //int outputs[numoutputs] = {ThicknessEnum};
    31                 int outputs[numoutputs] = {SurfaceEnum,VxEnum,VyEnum,VelEnum,BalancethicknessDiffusionCoefficientEnum};
     29                const int numoutputs = 4;
     30                int outputs[numoutputs] = {SurfaceEnum,VxEnum,VyEnum,VelEnum};
    3231                femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],numoutputs);
    3332        }

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

    r18476 r18565  
    2929        iomodel->FetchDataToInput(elements,InversionVyObsEnum,0.);
    3030        iomodel->FetchDataToInput(elements,InversionThicknessObsEnum,0.);
     31        iomodel->FetchDataToInput(elements,InversionSurfaceObsEnum,0.);
    3132
    3233        iomodel->FetchData(5,InversionControlParametersEnum,InversionCostFunctionsEnum,InversionCostFunctionsCoefficientsEnum,InversionMinParametersEnum,InversionMaxParametersEnum);
    … …  
    4243                        case FrictionCoefficientEnum:
    4344                        case BalancethicknessApparentMassbalanceEnum:
     45                        case BalancethicknessOmegaEnum:
    4446                                iomodel->FetchData(1,control);
    4547                                break;

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r18553 r18565  
    165165        InversionStepThresholdEnum,
    166166        InversionThicknessObsEnum,
     167        InversionSurfaceObsEnum,
    167168        InversionVxObsEnum,
    168169        InversionVyObsEnum,
    … …  
    307308        BalancethicknessDiffusionCoefficientEnum,
    308309        BalancethicknessCmuEnum,
     310        BalancethicknessOmegaEnum,
     311        BalancethicknessD0Enum,
    309312        /*}}}*/
    310313        /*Surfaceforcings{{{*/
    … …  
    527530        TemperaturePicardEnum,
    528531        ThicknessAbsMisfitEnum,
     532        SurfaceAbsMisfitEnum,
    529533        VelEnum,
    530534        VelocityEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r18553 r18565  
    173173                case InversionStepThresholdEnum : return "InversionStepThreshold";
    174174                case InversionThicknessObsEnum : return "InversionThicknessObs";
     175                case InversionSurfaceObsEnum : return "InversionSurfaceObs";
    175176                case InversionVxObsEnum : return "InversionVxObs";
    176177                case InversionVyObsEnum : return "InversionVyObs";
    … …  
    315316                case BalancethicknessDiffusionCoefficientEnum : return "BalancethicknessDiffusionCoefficient";
    316317                case BalancethicknessCmuEnum : return "BalancethicknessCmu";
     318                case BalancethicknessOmegaEnum : return "BalancethicknessOmega";
     319                case BalancethicknessD0Enum : return "BalancethicknessD0";
    317320                case SurfaceforcingsEnum : return "Surfaceforcings";
    318321                case SMBEnum : return "SMB";
    … …  
    518521                case TemperaturePicardEnum : return "TemperaturePicard";
    519522                case ThicknessAbsMisfitEnum : return "ThicknessAbsMisfit";
     523                case SurfaceAbsMisfitEnum : return "SurfaceAbsMisfit";
    520524                case VelEnum : return "Vel";
    521525                case VelocityEnum : return "Velocity";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r18553 r18565  
    176176              else if (strcmp(name,"InversionStepThreshold")==0) return InversionStepThresholdEnum;
    177177              else if (strcmp(name,"InversionThicknessObs")==0) return InversionThicknessObsEnum;
     178              else if (strcmp(name,"InversionSurfaceObs")==0) return InversionSurfaceObsEnum;
    178179              else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum;
    179180              else if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum;
    … …  
    259260              else if (strcmp(name,"QmuIsdakota")==0) return QmuIsdakotaEnum;
    260261              else if (strcmp(name,"MassFluxSegments")==0) return MassFluxSegmentsEnum;
    261               else if (strcmp(name,"MassFluxSegmentsPresent")==0) return MassFluxSegmentsPresentEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"QmuMassFluxSegmentsPresent")==0) return QmuMassFluxSegmentsPresentEnum;
     265              if (strcmp(name,"MassFluxSegmentsPresent")==0) return MassFluxSegmentsPresentEnum;
     266              else if (strcmp(name,"QmuMassFluxSegmentsPresent")==0) return QmuMassFluxSegmentsPresentEnum;
    266267              else if (strcmp(name,"QmuNumberofpartitions")==0) return QmuNumberofpartitionsEnum;
    267268              else if (strcmp(name,"QmuNumberofresponses")==0) return QmuNumberofresponsesEnum;
    … …  
    321322              else if (strcmp(name,"BalancethicknessDiffusionCoefficient")==0) return BalancethicknessDiffusionCoefficientEnum;
    322323              else if (strcmp(name,"BalancethicknessCmu")==0) return BalancethicknessCmuEnum;
     324              else if (strcmp(name,"BalancethicknessOmega")==0) return BalancethicknessOmegaEnum;
     325              else if (strcmp(name,"BalancethicknessD0")==0) return BalancethicknessD0Enum;
    323326              else if (strcmp(name,"Surfaceforcings")==0) return SurfaceforcingsEnum;
    324327              else if (strcmp(name,"SMB")==0) return SMBEnum;
    … …  
    380383              else if (strcmp(name,"HydrologyShreveAnalysis")==0) return HydrologyShreveAnalysisEnum;
    381384              else if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum;
    382               else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum;
    383               else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum;
    384               else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum;
     388              if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum;
     389              else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum;
     390              else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
     391              else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum;
    389392              else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
    390393              else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum;
    … …  
    503506              else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
    504507              else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
    505               else if (strcmp(name,"Converged")==0) return ConvergedEnum;
    506               else if (strcmp(name,"Fill")==0) return FillEnum;
    507               else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"Friction")==0) return FrictionEnum;
     511              if (strcmp(name,"Converged")==0) return ConvergedEnum;
     512              else if (strcmp(name,"Fill")==0) return FillEnum;
     513              else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum;
     514              else if (strcmp(name,"Friction")==0) return FrictionEnum;
    512515              else if (strcmp(name,"Internal")==0) return InternalEnum;
    513516              else if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
    … …  
    530533              else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
    531534              else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum;
     535              else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum;
    532536              else if (strcmp(name,"Vel")==0) return VelEnum;
    533537              else if (strcmp(name,"Velocity")==0) return VelocityEnum;
    … …  
    625629              else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum;
    626630              else if (strcmp(name,"MassfluxatgateName")==0) return MassfluxatgateNameEnum;
    627               else if (strcmp(name,"MassfluxatgateSegments")==0) return MassfluxatgateSegmentsEnum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"MassfluxatgateSegments")==0) return MassfluxatgateSegmentsEnum;
    628635              else if (strcmp(name,"MisfitName")==0) return MisfitNameEnum;
    629636              else if (strcmp(name,"MisfitModelEnum")==0) return MisfitModelEnumEnum;
    630637              else if (strcmp(name,"MisfitObservation")==0) return MisfitObservationEnum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"MisfitObservationEnum")==0) return MisfitObservationEnumEnum;
     638              else if (strcmp(name,"MisfitObservationEnum")==0) return MisfitObservationEnumEnum;
    635639              else if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum;
    636640              else if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum;
    … …  
    748752              else if (strcmp(name,"SurfaceforcingsAirQuadDragCoef")==0) return SurfaceforcingsAirQuadDragCoefEnum;
    749753              else if (strcmp(name,"SurfaceforcingsWindVx")==0) return SurfaceforcingsWindVxEnum;
    750               else if (strcmp(name,"SurfaceforcingsWindVy")==0) return SurfaceforcingsWindVyEnum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"SurfaceforcingsWindVy")==0) return SurfaceforcingsWindVyEnum;
    751758              else if (strcmp(name,"Matseaice")==0) return MatseaiceEnum;
    752759              else if (strcmp(name,"MaterialsPoisson")==0) return MaterialsPoissonEnum;
    753760              else if (strcmp(name,"MaterialsYoungModulus")==0) return MaterialsYoungModulusEnum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"MaterialsTimeRelaxationStress")==0) return MaterialsTimeRelaxationStressEnum;
     761              else if (strcmp(name,"MaterialsTimeRelaxationStress")==0) return MaterialsTimeRelaxationStressEnum;
    758762              else if (strcmp(name,"MaterialsTimeRelaxationDamage")==0) return MaterialsTimeRelaxationDamageEnum;
    759763              else if (strcmp(name,"MaterialsRidgingExponent")==0) return MaterialsRidgingExponentEnum;

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

    r18555 r18565  
    1010                stabilization     = 0;
    1111
    12                 Cmu               = NaN;
     12                omega             = NaN;
    1313        end
    1414        methods
    … …  
    3535                        md = checkfield(md,'fieldname','balancethickness.stabilization','size',[1 1],'values',[0 1 2 3]);
    3636
    37                         %md = checkfield(md,'fieldname','balancethickness.Cmu','size',[md.mesh.numberofvertices 1],'NaN',1,'>=',0);
     37                        %md = checkfield(md,'fieldname','balancethickness.omega','size',[md.mesh.numberofvertices 1],'NaN',1,'>=',0);
    3838                end % }}}
    3939                function disp(obj) % {{{
    … …  
    5353                        WriteData(fid,'object',obj,'fieldname','stabilization','format','Integer');
    5454
    55                         WriteData(fid,'object',obj,'fieldname','Cmu','format','DoubleMat','mattype',1);
     55                        WriteData(fid,'object',obj,'fieldname','omega','format','DoubleMat','mattype',1);
    5656                end % }}}
    5757        end

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

    r18476 r18565  
    2222                vel_obs                     = NaN
    2323                thickness_obs               = NaN
     24                surface_obs               = NaN
    2425
    2526        end
    … …  
    7273                        md = checkfield(md,'fieldname','inversion.control_parameters','cell',1,'values',...
    7374                                {'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'DamageDbar',...
    74                                 'Vx' 'Vy' 'Thickness' 'BalancethicknessNu' 'BalancethicknessApparentMassbalance'});
     75                                'Vx' 'Vy' 'Thickness' 'BalancethicknessOmega' 'BalancethicknessApparentMassbalance'});
    7576                        md = checkfield(md,'fieldname','inversion.maxsteps','numel',1,'>=',0);
    7677                        md = checkfield(md,'fieldname','inversion.maxiter','numel',1,'>=',0);
    … …  
    8485                        if solution==BalancethicknessSolutionEnum()
    8586                                md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
     87                                md = checkfield(md,'fieldname','inversion.surface_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
    8688                        elseif solution==BalancethicknessSoftSolutionEnum()
    8789                                md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
    … …  
    110112                        fielddisplay(obj,'vel_obs','observed velocity magnitude [m/yr]');
    111113                        fielddisplay(obj,'thickness_obs','observed thickness [m]');
     114                        fielddisplay(obj,'surface_obs','observed surface elevation [m]');
    112115                        disp('Available cost functions:');
    113116                        disp('   101: SurfaceAbsVelMisfit');
    … …  
    145148                        end
    146149                        WriteData(fid,'object',obj,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',mattype);
     150                        WriteData(fid,'object',obj,'class','inversion','fieldname','surface_obs','format','DoubleMat','mattype',mattype);
    147151
    148152                        %process control parameters
    … …  
    170174                        pos=find(obj.cost_functions==505); data(pos)=ThicknessAcrossGradientEnum();
    171175                        pos=find(obj.cost_functions==506); data(pos)=BalancethicknessMisfitEnum();
     176                        pos=find(obj.cost_functions==601); data(pos)=SurfaceAbsMisfitEnum();
    172177                        WriteData(fid,'data',data,'enum',InversionCostFunctionsEnum(),'format','DoubleMat','mattype',3);
    173178                        WriteData(fid,'data',num_cost_functions,'enum',InversionNumCostFunctionsEnum(),'format','Integer');

  • issm/trunk-jpl/src/m/enum/EnumDefinitions.py

    r18553 r18565  
    165165def InversionStepThresholdEnum(): return StringToEnum("InversionStepThreshold")[0]
    166166def InversionThicknessObsEnum(): return StringToEnum("InversionThicknessObs")[0]
     167def InversionSurfaceObsEnum(): return StringToEnum("InversionSurfaceObs")[0]
    167168def InversionVxObsEnum(): return StringToEnum("InversionVxObs")[0]
    168169def InversionVyObsEnum(): return StringToEnum("InversionVyObs")[0]
    … …  
    307308def BalancethicknessDiffusionCoefficientEnum(): return StringToEnum("BalancethicknessDiffusionCoefficient")[0]
    308309def BalancethicknessCmuEnum(): return StringToEnum("BalancethicknessCmu")[0]
     310def BalancethicknessOmegaEnum(): return StringToEnum("BalancethicknessOmega")[0]
     311def BalancethicknessD0Enum(): return StringToEnum("BalancethicknessD0")[0]
    309312def SurfaceforcingsEnum(): return StringToEnum("Surfaceforcings")[0]
    310313def SMBEnum(): return StringToEnum("SMB")[0]
    … …  
    510513def TemperaturePicardEnum(): return StringToEnum("TemperaturePicard")[0]
    511514def ThicknessAbsMisfitEnum(): return StringToEnum("ThicknessAbsMisfit")[0]
     515def SurfaceAbsMisfitEnum(): return StringToEnum("SurfaceAbsMisfit")[0]
    512516def VelEnum(): return StringToEnum("Vel")[0]
    513517def VelocityEnum(): return StringToEnum("Velocity")[0]
Note: See TracChangeset for help on using the changeset viewer.