![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
- Timestamp:
- 07/10/20 17:15:11 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/modules/InputUpdateFromDakotax/InputUpdateFromDakotax.cpp
r25233 r25258 13 13 #include "../InputUpdateFromVectorDakotax/InputUpdateFromVectorDakotax.h" 14 14 15 void InputUpdateFromDakotax(FemModel* femmodel,double* variables,char* *variables_descriptors,int numdakotavariables){ /*{{{*/15 void InputUpdateFromDakotax(FemModel* femmodel,double* variables,char* *variables_descriptors,int numdakotavariables){ /*{{{*/ 16 16 17 17 int i,j,k,l; … … 46 46 47 47 /*Go through all dakota descriptors, ex: "rho_ice","thermal_conductivity","thickness1","thickness2", etc ..., and 48 * for each descriptor, take the variable value and plug it into the inputs: */ 48 * for each descriptor, take the variable value and plug it into the inputs (more or less :)): 49 * We also start with distributed and standard values , as they tend to be used to pluck data from a multi-modle ensemble (mme) 50 * which can then be scaled. Doing the scaling first would be impractical, as the entire mme would have to be scaled, 51 * which is a waste of time:*/ 52 53 variablecount=0; 49 54 for(i=0;i<numdakotavariables;i++){ //these are the dakota variables, for all partitions. 50 55 51 56 descriptor=variables_descriptors[i]; 52 57 53 if (VerboseQmu())_printf0_(" updating variable " << descriptor << "\n"); 54 55 /*From descriptor, figure out if the variable is scaled, indexed, nodal, or just a simple variable: */ 58 59 /*From descriptor, figure out if the variable is scaled, indexed, distributed or just a simple variable: */ 60 if (strncmp(descriptor,"scaled_",7)==0){ 61 /*we are skipping these for now.*/ 62 variable_partition=variable_partitions[variablecount]; 63 npart=variable_partitions_npart[variablecount]; 64 nt=variable_partitions_nt[variablecount]; 65 66 /*increment i to skip the distributed values just collected: */ 67 i+=npart*nt-1; //careful, the for loop will add 1. 68 } 69 else if (strncmp(descriptor,"indexed_",8)==0){ 70 /*we are skipping these for now.*/ 71 } 72 else if (strncmp(descriptor,"nodal_",8)==0){ 73 /*we are skipping these for now.*/ 74 } 75 76 else if (strncmp(descriptor,"distributed_",12)==0){ 77 78 if (VerboseQmu())_printf0_(" updating variable " << descriptor << "\n"); 79 80 /*recover partition vector: */ 81 variable_partition=variable_partitions[variablecount]; 82 npart=variable_partitions_npart[variablecount]; 83 84 /*Variable is distributed. Determine root name of variable (ex: distributed_DragCoefficient_1 -> DragCoefficient). 85 * Allocate distributed_values and fill the distributed_values with the next npart variables: */ 86 87 memcpy(root,strstr(descriptor,"_")+1,(strlen(strstr(descriptor,"_")+1)+1)*sizeof(char)); 88 *strstr(root,"_")='\0'; 89 90 distributed_values=xNew<double>(npart); 91 for(j=0;j<npart;j++){ 92 distributed_values[j]=variables[i+j]; 93 } 94 95 //for (int j=0;j<npart;j++)_printf_(j << ":" << distributed_values[j] << "\n"); 96 97 //Call specialty code: 98 InputUpdateSpecialtyCode(femmodel,distributed_values,variable_partition,npart,root); 99 100 /*increment i to skip the distributed values just collected: */ 101 i+=npart-1; //careful, the for loop will add 1. 102 103 /*Free allocations: */ 104 xDelete<double>(parameter); 105 xDelete<double>(distributed_values); 106 } 107 else{ 108 /*Ok, standard variable, just update inputs using the variable: */ 109 if (VerboseQmu())_printf0_(" updating variable " << descriptor << "\n"); 110 InputUpdateFromConstantx(femmodel,variables[i],StringToEnumx(descriptor)); 111 } 112 variablecount++; 113 } 114 115 variablecount=0; 116 /*now deal with scaled variabes:*/ 117 for(i=0;i<numdakotavariables;i++){ //these are the dakota variables, for all partitions. 118 119 descriptor=variables_descriptors[i]; 120 121 122 /*From descriptor, figure out if the variable is scaled, indexed, distributed or just a simple variable: */ 56 123 if (strncmp(descriptor,"scaled_",7)==0){ 124 125 if (VerboseQmu())_printf0_(" updating variable " << descriptor << "\n"); 57 126 58 127 /*recover partition vector: */ … … 61 130 nt=variable_partitions_nt[variablecount]; 62 131 63 /*Variable is scaled. Determine root name of variable (ex: scaled_DragCoefficient_1 -> DragCoefficient). Allocate distributed_values and fill the 64 * distributed_values with the next npart variables: */ 65 66 //strcpy(root,strstr(descriptor,"_")+1); *strstr(root,"_")='\0'; 132 /* Variable is scaled, determine its root name (ex: scaled_DragCoefficient_1 -> DragCoefficient). Allocate distributed_values and fill the 133 * distributed_values with the next npart variables coming from Dakota: */ 67 134 memcpy(root,strstr(descriptor,"_")+1,(strlen(strstr(descriptor,"_")+1)+1)*sizeof(char)); 68 135 *strstr(root,"_")='\0'; … … 74 141 75 142 76 /*Now, pick up the parameter corresponding to root: */ 77 femmodel->parameters->FindParamInDataset(¶meter,&nrows,&ncols,QmuVariableDescriptorsEnum,StringToEnumx(root)); 78 79 /*We've got the parameter, we need to update it using the partition vector, and the distributed_values. 80 In addition, the parameter can be either a static or transient (nrows+1) vector. Finally, the partition vectors can include 81 -1 (meaning, don't update). */ 82 83 //_printf_("nrows: " << nrows << " numberofvertices: " << numberofvertices << " numberofelements: " << numberofelements << "\n"); 84 85 if(ncols!=nt){ 86 /*we are trying to update a col sized transient input by scaling with a matrix of col size nt. This can only work if nt==1, otherwise, error out: */ 87 if (nt!=1) _error_("InputUpdateFromDakotax error message: transient input being updated should be the same col size as the number of time step in the qmu variable specificationi"); 88 } 89 90 if(nt==1){ 91 /*scale all the columns by the same vector:*/ 92 if (nrows==numberofvertices || nrows==(numberofvertices+1)){ 93 for(k=0;k<numberofvertices;k++){ 94 if (variable_partition[k]==-1)continue; 95 else{ 96 for(l=0;l<ncols;l++){ 97 *(parameter+ncols*k+l)=*(parameter+ncols*k+l)*distributed_values[(int)variable_partition[k]]; 143 int inputscaling=0; 144 if(!inputscaling){ 145 /*Scale variable outside of the inputs:{{{*/ 146 147 /*Now, pick up the parameter corresponding to root: */ 148 femmodel->parameters->FindParamInDataset(¶meter,&nrows,&ncols,QmuVariableDescriptorsEnum,StringToEnumx(root)); 149 150 /*We've got the parameter, we need to update it using the partition vector, and the distributed_values. 151 In addition, the parameter can be either a static or transient (nrows+1) vector. Finally, the partition vectors can include 152 -1 (meaning, don't update). */ 153 154 //_printf_("nrows: " << nrows << " numberofvertices: " << numberofvertices << " numberofelements: " << numberofelements << "\n"); 155 156 if(ncols!=nt){ 157 /*we are trying to update a col sized transient input by scaling with a matrix of col size nt. This can only work if nt==1, otherwise, error out: */ 158 if (nt!=1) _error_("InputUpdateFromDakotax error message: transient input being updated should be the same col size as the number of time step in the qmu variable specificationi"); 159 } 160 161 if(nt==1){ 162 /*scale all the columns by the same vector:*/ 163 if (nrows==numberofvertices || nrows==(numberofvertices+1)){ 164 for(k=0;k<numberofvertices;k++){ 165 if (variable_partition[k]==-1)continue; 166 else{ 167 for(l=0;l<ncols;l++){ 168 *(parameter+ncols*k+l)=*(parameter+ncols*k+l)*distributed_values[(int)variable_partition[k]]; 169 } 98 170 } 99 171 } 100 172 } 173 else if (nrows==numberofelements || nrows==(numberofelements+1)){ 174 for(k=0;k<numberofelements;k++){ 175 if (variable_partition[k]==-1)continue; 176 else{ 177 for(l=0;l<ncols;l++){ 178 *(parameter+ncols*k+l)=*(parameter+ncols*k+l)*distributed_values[(int)variable_partition[k]]; 179 } 180 } 181 } 182 183 } 184 else _error_("partitioning vector should be either elements or vertex sized!"); 185 101 186 } 102 else if (nrows==numberofelements || nrows==(numberofelements+1)){ 103 for(k=0;k<numberofelements;k++){ 104 if (variable_partition[k]==-1)continue; 105 else{ 106 for(l=0;l<ncols;l++){ 107 *(parameter+ncols*k+l)=*(parameter+ncols*k+l)*distributed_values[(int)variable_partition[k]]; 187 else{ 188 /*scale all the columns by the scalar matrix:*/ 189 if (nrows==numberofvertices || nrows==(numberofvertices+1)){ 190 for(k=0;k<numberofvertices;k++){ 191 if (variable_partition[k]==-1)continue; 192 else{ 193 for(l=0;l<ncols;l++){ 194 *(parameter+ncols*k+l)=*(parameter+ncols*k+l)*distributed_values[(int)variable_partition[k]*nt+l]; 195 } 108 196 } 109 197 } 110 198 } 111 112 } 113 else _error_("partitioning vector should be either elements or vertex sized!"); 114 115 } 116 else{ 117 /*scale all the columns by the scalar matrix:*/ 118 if (nrows==numberofvertices || nrows==(numberofvertices+1)){ 119 for(k=0;k<numberofvertices;k++){ 120 if (variable_partition[k]==-1)continue; 121 else{ 122 for(l=0;l<ncols;l++){ 123 *(parameter+ncols*k+l)=*(parameter+ncols*k+l)*distributed_values[(int)variable_partition[k]*nt+l]; 199 else if (nrows==numberofelements || nrows==(numberofelements+1)){ 200 for(k=0;k<numberofelements;k++){ 201 if (variable_partition[k]==-1)continue; 202 else{ 203 for(l=0;l<ncols;l++){ 204 *(parameter+ncols*k+l)=*(parameter+ncols*k+l)*distributed_values[(int)variable_partition[k]*nt+l]; 205 } 124 206 } 125 207 } 208 126 209 } 210 else _error_("partitioning vector should be either elements or vertex sized!"); 127 211 } 128 else if (nrows==numberofelements || nrows==(numberofelements+1)){ 129 for(k=0;k<numberofelements;k++){ 130 if (variable_partition[k]==-1)continue; 131 else{ 132 for(l=0;l<ncols;l++){ 133 *(parameter+ncols*k+l)=*(parameter+ncols*k+l)*distributed_values[(int)variable_partition[k]*nt+l]; 134 } 135 } 136 } 137 138 } 139 else _error_("partitioning vector should be either elements or vertex sized!"); 140 } 141 142 #ifdef _DEBUG_ 143 PetscSynchronizedPrintf(IssmComm::GetComm(),"Parameter matrix:"); 144 PetscSynchronizedFlush(IssmComm::GetComm()); 145 for(l=0;l<ncols;l++){ 146 PetscSynchronizedPrintf(IssmComm::GetComm()," time %i\n",l); 147 PetscSynchronizedFlush(IssmComm::GetComm()); 148 149 for(k=0;k<numberofvertices;k++){ 150 PetscSynchronizedPrintf(IssmComm::GetComm()," node %i value %g\n",k+1,*(parameter+k*ncols+l)); 151 PetscSynchronizedFlush(IssmComm::GetComm()); 152 } 153 } 154 PetscSynchronizedPrintf(IssmComm::GetComm()," descriptor: %s root %s enum: %i\n",descriptor,root,StringToEnumx(root)); 155 PetscSynchronizedFlush(IssmComm::GetComm()); 156 #endif 157 158 /*Update inputs using the parameter matrix: */ 159 if(nrows==numberofvertices || (nrows==numberofvertices+1)) 160 InputUpdateFromMatrixDakotax(femmodel, parameter, nrows,ncols,StringToEnumx(root), VertexEnum); 161 else 162 InputUpdateFromMatrixDakotax(femmodel, parameter, nrows,ncols,StringToEnumx(root), ElementEnum); 212 213 /*Update inputs using the parameter matrix: */ 214 if(nrows==numberofvertices || (nrows==numberofvertices+1)) 215 InputUpdateFromMatrixDakotax(femmodel, parameter, nrows,ncols,StringToEnumx(root), VertexEnum); 216 else 217 InputUpdateFromMatrixDakotax(femmodel, parameter, nrows,ncols,StringToEnumx(root), ElementEnum); /*}}}*/ 218 } 219 else{ 220 /*Scale variable inside the inputs:*/ 221 InputScaleFromDakotax(femmodel, distributed_values, variable_partition,npart, nt, StringToEnumx(root)); 222 } 163 223 164 224 /*increment i to skip the distributed values just collected: */ … … 168 228 xDelete<double>(parameter); 169 229 xDelete<double>(distributed_values); 170 }171 else if (strncmp(descriptor,"distributed_",12)==0){172 173 /*recover partition vector: */174 variable_partition=variable_partitions[variablecount];175 npart=variable_partitions_npart[variablecount];176 177 /*Variable is distributed. Determine root name of variable (ex: distributed_DragCoefficient_1 -> DragCoefficient).178 * Allocate distributed_values and fill the distributed_values with the next npart variables: */179 180 memcpy(root,strstr(descriptor,"_")+1,(strlen(strstr(descriptor,"_")+1)+1)*sizeof(char));181 *strstr(root,"_")='\0';182 183 distributed_values=xNew<double>(npart);184 for(j=0;j<npart;j++){185 distributed_values[j]=variables[i+j];186 }187 188 //for (int j=0;j<npart;j++)_printf_(j << ":" << distributed_values[j] << "\n");189 190 //Call specialty code:191 InputUpdateSpecialtyCode(femmodel,distributed_values,variable_partition,npart,root);192 193 /*increment i to skip the distributed values just collected: */194 i+=npart-1; //careful, the for loop will add 1.195 196 /*Free allocations: */197 xDelete<double>(parameter);198 xDelete<double>(distributed_values);199 }200 else if (strncmp(descriptor,"indexed_",8)==0){201 _error_("indexed variables not supported yet!");202 }203 else if (strncmp(descriptor,"nodal_",8)==0){204 _error_("nodal variables not supported yet!");205 }206 else{207 /*Ok, standard variable, just update inputs using the variable: */208 InputUpdateFromConstantx(femmodel,variables[i],StringToEnumx(descriptor));209 230 } 210 231 variablecount++; … … 230 251 231 252 if(VerboseQmu()){ 232 _printf0_(" UpdatingSurfaceloadModelid MME, with ids: ");253 _printf0_(" SurfaceloadModelid MME, with ids: "); 233 254 for (int i=0;i<npart;i++)_printf0_((int)distributed_values[i]+1 << " "); 234 255 _printf0_("\n"); … … 309 330 transientinput2->Configure(femmodel->parameters); 310 331 } //}}} 332 void InputScaleFromDakotax(FemModel* femmodel,IssmDouble* distributed_values,IssmDouble* partition, int npart, int nt, int name){ /*{{{*/ 333 334 /*Copy input:*/ 335 femmodel->inputs2->DuplicateInput(name,DummyEnum); 336 337 /*Go through elements, copy input name to dummy, and scale it using the distributed_values and the partition vector:*/ 338 for(int i=0;i<femmodel->elements->Size();i++){ 339 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i)); 340 element->InputScaleFromDakota(distributed_values,partition,npart,nt,name); 341 } 342 343 /*We created a dummy input, which was a scaled copy of the name input. Now wipe 344 * out the name input with the new input:*/ 345 femmodel->inputs2->ChangeEnum(DummyEnum,name); 346 347 348 } /*}}}*/
Note:
See TracChangeset
for help on using the changeset viewer.
