Changeset 25273


Ignore:
Timestamp:
07/13/20 22:04:50 (5 years ago)
Author:
Eric.Larour
Message:

CHG: dslmodelid should be a double param for uq runs.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r25219 r25273  
    135135                xDelete<int>(pN);
    136136                /*}}}*/
    137                 /*now do the same with the dsl.sea_surface_height_change_above_geoid field:{{{ */
    138                 iomodel->FetchData(&pstr,&pM,&pN,&nummodels,"md.dsl.sea_surface_height_change_above_geoid");
    139 
    140                 for (int i=0;i<nummodels;i++){
    141                         M=pM[i];
    142                         N=pN[i];
    143                         str=pstr[i];
    144                
    145 
    146                         //recover time vector:
    147                         times=xNew<IssmDouble>(N);
    148                         for(int t=0;t<N;t++) times[t] = str[(M-1)*N+t];
    149                        
    150                         TransientInput2* transientinput=inputs2->SetDatasetTransientInput(DslSeaSurfaceHeightChangeAboveGeoidEnum,i, times,N);
    151        
    152                         for(int j=0;j<elements->Size();j++){
    153 
    154                                 /*Get the right transient input*/
    155                                 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
    156 
    157                                 /*Get values and lid list*/
    158                                 const int   numvertices = element->GetNumberOfVertices();
    159                                 int        *vertexlids = xNew<int>(numvertices);
    160                                 int        *vertexsids = xNew<int>(numvertices);
    161 
    162 
    163                                 /*Recover vertices ids needed to initialize inputs*/
    164                                 _assert_(iomodel->elements);
    165                                 for(int k=0;k<numvertices;k++){
    166                                         vertexsids[k] =reCast<int>(iomodel->elements[numvertices*element->Sid()+k]-1); //ids for vertices are in the elements array from Matlab
    167                                         vertexlids[k]=iomodel->my_vertices_lids[vertexsids[k]];
    168                                 }
    169 
    170                                 IssmDouble* values=xNew<IssmDouble>(numvertices);
    171 
    172                                 for(int t=0;t<N;t++){
    173                                         for (int k=0;k<numvertices;k++)values[k]=str[N*vertexsids[k]+t];
    174 
    175                                         switch(element->ObjectEnum()){
    176                                                 case TriaEnum:  transientinput->AddTriaTimeInput( t,numvertices,vertexlids,values,P1Enum); break;
    177                                                 case PentaEnum: transientinput->AddPentaTimeInput(t,numvertices,vertexlids,values,P1Enum); break;
    178                                                 default: _error_("Not implemented yet");
    179                                         }
    180                                 }
    181                                 xDelete<IssmDouble>(values);
    182                                 xDelete<int>(vertexlids);
    183                                 xDelete<int>(vertexsids);
    184                         }
    185                        
    186                         xDelete<IssmDouble>(times);
    187                 }
    188                
    189                 /*Delete data:*/
    190                 for(int i=0;i<nummodels;i++){
    191                         IssmDouble* str=pstr[i];
    192                         xDelete<IssmDouble>(str);
    193                 }
    194                 xDelete<IssmDouble*>(pstr);
    195                 xDelete<int>(pM);
    196                 xDelete<int>(pN);
    197                 /*}}}*/
    198                 /*now do the same with the dsl.sea_water_pressure_change_at_sea_floor field:{{{ */
    199                 iomodel->FetchData(&pstr,&pM,&pN,&nummodels,"md.dsl.sea_water_pressure_change_at_sea_floor");
    200 
    201                 for (int i=0;i<nummodels;i++){
    202                         M=pM[i];
    203                         N=pN[i];
    204                         str=pstr[i];
    205 
    206                         //recover time vector:
    207                         times=xNew<IssmDouble>(N);
    208                         for(int t=0;t<N;t++) times[t] = str[(M-1)*N+t];
    209 
    210                         TransientInput2* transientinput=inputs2->SetDatasetTransientInput(DslSeaWaterPressureChangeAtSeaFloor,i, times,N);
    211        
    212                         for(int j=0;j<elements->Size();j++){
    213 
    214                                 /*Get the right transient input*/
    215                                 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
    216 
    217                                 /*Get values and lid list*/
    218                                 const int   numvertices = element->GetNumberOfVertices();
    219                                 int        *vertexlids = xNew<int>(numvertices);
    220                                 int        *vertexsids = xNew<int>(numvertices);
    221                                
    222                                
    223                                 /*Recover vertices ids needed to initialize inputs*/
    224                                 _assert_(iomodel->elements);
    225                                 for(int k=0;k<numvertices;k++){
    226                                         vertexsids[k] =reCast<int>(iomodel->elements[numvertices*element->Sid()+k]-1); //ids for vertices are in the elements array from Matlab
    227                                         vertexlids[k]=iomodel->my_vertices_lids[vertexsids[k]];
    228                                 }
    229                                
    230                                 IssmDouble* values=xNew<IssmDouble>(numvertices);
    231 
    232                                 for(int t=0;t<N;t++){
    233                                         for (int k=0;k<numvertices;k++)values[k]=str[N*vertexsids[k]+t];
    234 
    235                                         switch(element->ObjectEnum()){
    236                                                 case TriaEnum:  transientinput->AddTriaTimeInput( t,numvertices,vertexlids,values,P1Enum); break;
    237                                                 case PentaEnum: transientinput->AddPentaTimeInput(t,numvertices,vertexlids,values,P1Enum); break;
    238                                                 default: _error_("Not implemented yet");
    239                                         }
    240                                 }
    241                                 xDelete<IssmDouble>(values);
    242                                 xDelete<int>(vertexlids);
    243                                 xDelete<int>(vertexsids);
    244                         }
    245                         xDelete<IssmDouble>(times);
    246                 }
    247                
    248                 /*Delete data:*/
    249                 for(int i=0;i<nummodels;i++){
    250                         IssmDouble* str=pstr[i];
    251                         xDelete<IssmDouble>(str);
    252                 }
    253                 xDelete<IssmDouble*>(pstr);
    254                 xDelete<int>(pM);
    255                 xDelete<int>(pN);
    256                 /*}}}*/
     137                iomodel->FetchDataToDatasetInput(inputs2,elements,"md.dsl.sea_surface_height_change_above_geoid",DslSeaSurfaceHeightChangeAboveGeoidEnum);
     138                iomodel->FetchDataToDatasetInput(inputs2,elements,"md.dsl.sea_water_pressure_change_at_sea_floor",DslSeaWaterPressureChangeAtSeaFloor);
    257139
    258140        } /*}}}*/
    … …  
    329211        parameters->AddObject(iomodel->CopyConstantObject("md.dsl.compute_fingerprints",DslComputeFingerprintsEnum));
    330212        if(dslmodel==2){
    331                 int modelid;
     213                IssmDouble modelid;
    332214                int nummodels;
    333 
    334                 parameters->AddObject(iomodel->CopyConstantObject("md.dsl.modelid",DslModelidEnum));
     215               
     216                /*create double param, not int param, because Dakota will be updating it as
     217                 * a double potentially: */
     218                iomodel->FetchData(&modelid,"md.dsl.modelid");
     219                parameters->AddObject(new DoubleParam(DslModelEnum,modelid));
    335220                parameters->AddObject(iomodel->CopyConstantObject("md.dsl.nummodels",DslNummodelsEnum));
    336                 iomodel->FetchData(&modelid,"md.dsl.modelid");
    337221                iomodel->FetchData(&nummodels,"md.dsl.nummodels");
    338222
Note: See TracChangeset for help on using the changeset viewer.