![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 13432
- Timestamp:
- 09/24/12 22:26:49 (12 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 7 edited
-
Makefile.am (modified) (2 diffs)
-
classes/IoModel.cpp (modified) (27 diffs)
-
classes/IoModel.h (modified) (1 diff)
-
modules/AutodiffDriversx/AutodiffDriversx.cpp (modified) (5 diffs)
-
modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp (modified) (2 diffs)
-
modules/ModelProcessorx/ModelProcessorx.cpp (modified) (2 diffs)
-
modules/RequestedDependentsx/RequestedDependentsx.cpp (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/Makefile.am
r13325 r13432 54 54 ./classes/objects/Constraints/SpcTransient.cpp\ 55 55 ./classes/objects/Constraints/SpcTransient.h\ 56 ./classes/objects/IndependentObject.h\ 57 ./classes/objects/IndependentObject.cpp\ 58 ./classes/objects/DependentObject.h\ 59 ./classes/objects/DependentObject.cpp\ 56 60 ./classes/DofIndexing.h\ 57 61 ./classes/DofIndexing.cpp\ … … 152 156 ./classes/objects/Params/TransientParam.h\ 153 157 ./classes/objects/Params/TransientParam.cpp\ 158 ./classes/objects/Params/DataSetParam.h\ 159 ./classes/objects/Params/DataSetParam.cpp\ 154 160 ./Container/Container.h\ 155 161 ./Container/Constraints.h\ -
issm/trunk-jpl/src/c/classes/IoModel.cpp
r13413 r13432 24 24 /*FUNCTION IoModel::IoModel(){{{*/ 25 25 IoModel::IoModel(){ 26 this->fid =NULL;27 this->data =NULL;28 this->independents =NULL;29 this->constants =NULL;30 31 this->my_elements =NULL;32 this->my_nodes =NULL;33 this->my_vertices =NULL;34 this->singlenodetoelementconnectivity =NULL;35 this->numbernodetoelementconnectivity =NULL;36 37 this->nodecounter =0;38 this->loadcounter =0;39 this->constraintcounter =0;26 this->fid=NULL; 27 this->data=NULL; 28 this->independents=NULL; 29 this->constants=NULL; 30 31 this->my_elements=NULL; 32 this->my_nodes=NULL; 33 this->my_vertices=NULL; 34 this->singlenodetoelementconnectivity=NULL; 35 this->numbernodetoelementconnectivity=NULL; 36 37 this->nodecounter=0; 38 this->loadcounter=0; 39 this->constraintcounter=0; 40 40 } 41 41 /*}}}*/ 42 42 /*FUNCTION IoModel::IoModel(FILE* iomodel_handle){{{*/ 43 43 IoModel::IoModel(FILE* iomodel_handle){ 44 44 45 45 /*First, keep track of the file handle: */ 46 46 this->fid=iomodel_handle; … … 48 48 /*Check that Enums are Synchronized*/ 49 49 this->CheckEnumSync(); 50 51 /*Initialize data: */ 52 this->data=xNew<IssmDouble*>(MaximumNumberOfEnums); 53 for(int i=0;i<MaximumNumberOfEnums;i++) this->data[i]=NULL; 54 55 /*If we are running in AD mode, we need to declare our independent variables now, 56 *and prevent them from being erased during successive calls to iomodel->FetchConstants, iomodel->FetchData and 57 iomodel->DeleteData:*/ 58 this->DeclareIndependents(); 50 59 51 60 /*Initialize and read constants:*/ … … 53 62 this->FetchConstants(); /*this routine goes through the input file, and fetches bools, ints, IssmDoubles and strings only, nothing memory intensive*/ 54 63 55 /*Initialize data: */56 this->data=xNew<IssmDouble*>(MaximumNumberOfEnums);57 for(int i=0;i<MaximumNumberOfEnums;i++) this->data[i]=NULL;58 59 /*Initialize array detecting whether data[i] is an independent AD mode variable: */60 this->independents=xNew<bool>(MaximumNumberOfEnums);61 for(int i=0;i<MaximumNumberOfEnums;i++) this->independents[i]=false;62 63 64 /*Initialize permanent data: */ 64 this->my_elements =NULL;65 this->my_nodes =NULL;66 this->my_vertices =NULL;67 this->singlenodetoelementconnectivity =NULL;68 this->numbernodetoelementconnectivity =NULL;69 70 this->nodecounter =0;71 this->loadcounter =0;72 this->constraintcounter =0;65 this->my_elements=NULL; 66 this->my_nodes=NULL; 67 this->my_vertices=NULL; 68 this->singlenodetoelementconnectivity=NULL; 69 this->numbernodetoelementconnectivity=NULL; 70 71 this->nodecounter=0; 72 this->loadcounter=0; 73 this->constraintcounter=0; 73 74 } 74 75 /*}}}*/ … … 91 92 xDelete<IssmDouble*>(this->data); 92 93 xDelete<bool>(this->independents); 94 delete this->independent_objects; 93 95 xDelete<bool>(this->my_elements); 94 96 xDelete<bool>(this->my_nodes); … … 197 199 } 198 200 /*}}}*/ 201 /*FUNCTION IoModel::DeclareIndependents{{{*/ 202 void IoModel::DeclareIndependents(void){ 203 204 205 int i; 206 bool autodiff = false; 207 int num_independent_objects; 208 209 int *names = NULL; 210 int *types = NULL; 211 212 int numberofvertices; 213 int dummy; 214 215 216 /*Initialize array detecting whether data[i] is an independent AD mode variable: */ 217 this->independents=xNew<bool>(MaximumNumberOfEnums); 218 for(i=0;i<MaximumNumberOfEnums;i++) this->independents[i]=false; 219 220 this->FetchData(&autodiff,AutodiffIsautodiffEnum); 221 if(autodiff){ 222 223 this->FetchData(&numberofvertices,MeshNumberofverticesEnum); 224 225 #ifdef _HAVE_ADOLC_ 226 /*build dataset made of independent objects: */ 227 this->independent_objects=new DataSet(); 228 this->FetchData(&num_independent_objects,AutodiffNumIndependentObjectsEnum); 229 if(num_independent_objects){ 230 this->FetchData(&names,&dummy,&dummy,AutodiffIndependentObjectNamesEnum); 231 this->FetchData(&types,&dummy,&dummy,AutodiffIndependentObjectTypesEnum); 232 233 /*create independent objects, and at the same time, fetch the corresponding independent variables, 234 *and declare them as such in ADOLC: */ 235 for(i=0;i<num_independent_objects;i++){ 236 237 IndependentObject* independent_object=NULL; 238 independent_object=new IndependentObject(names[i],types[i],numberofvertices); 239 240 /*add to independent_objects dataset:*/ 241 this->independent_objects->AddObject(independent_object); 242 243 /*now go fetch the independent variable: */ 244 independent_object->FetchIndependent(this); //supply the pointer to iomodel. 245 } 246 xDelete<int>(names); 247 xDelete<int>(types); 248 } 249 #else 250 /*if we asked for AD computations, we have a problem!: */ 251 _error_("Cannot carry out AD mode computations without support of ADOLC compiled in!"); 252 #endif 253 } 254 255 } 256 /*}}}*/ 199 257 /*FUNCTION IoModel::DeleteData(int num,...){{{*/ 200 258 void IoModel::DeleteData(int num,...){ 201 259 202 va_list ap; 203 int dataenum; 204 DoubleMatParam *parameter = NULL; 260 va_list ap; 261 int dataenum; 262 int i; 263 DoubleMatParam* parameter=NULL; 205 264 206 265 /*Go through the entire list of enums and delete the corresponding data from the iomodel-data dataset: */ 266 207 267 va_start(ap,num); 208 for(i nt i=0; i<num; i++){268 for(i = 0; i <num; i++){ 209 269 dataenum=va_arg(ap, int); 210 270 _assert_(dataenum<MaximumNumberOfEnums); … … 215 275 va_end(ap); 216 276 } /*}}}*/ 217 /*FUNCTION IoModel::DeleteData(IssmDouble* vector, int dataenum){{{*/277 /*FUNCTION IoModel::DeleteData(IssmDouble* {{{*/ 218 278 void IoModel::DeleteData(IssmDouble* vector, int dataenum){ 219 279 … … 226 286 227 287 extern int my_rank; 288 extern int num_procs; 228 289 229 290 /*record descriptions; */ 230 291 int record_enum; 231 292 int record_length; 232 int record_code; 293 int record_code; //1 to 7 number 233 294 234 295 /*records: */ 235 int booleanint = 0; 236 int integer = 0; 237 IssmPDouble scalar = 0; 238 char *string = NULL; 239 int string_size; 296 int booleanint=0; 297 int integer=0; 298 IssmPDouble pscalar=0; 299 IssmDouble scalar=0; 300 char* string=NULL; 301 int string_size; 240 302 241 303 /*Check that some fields have been allocated*/ … … 299 361 break; 300 362 case 3: 301 /*Read the scalar and broadcast it to other cpus:*/ 302 if(fread(&scalar,sizeof(IssmPDouble),1,this->fid)!=1) _error_("could not read scalar "); 303 #ifdef _HAVE_MPI_ 304 MPI_Bcast(&scalar,1,MPI_DOUBLE,0,MPI_COMM_WORLD); 305 #endif 363 /*Read the scalar and broadcast it to other cpus. However, if this record has already 364 * been read in DeclareIndependents, we should not re-read it, but grab it from the iomodel->data 365 * slots: */ 366 if(this->independents[record_enum]){ 367 scalar=*(this->data[record_enum]); 368 } 369 else{ 370 if(fread(&pscalar,sizeof(IssmPDouble),1,this->fid)!=1) _error_("could not read scalar "); 371 #ifdef _HAVE_MPI_ 372 MPI_Bcast(&pscalar,1,MPI_DOUBLE,0,MPI_COMM_WORLD); 373 #endif 374 scalar=reCast<IssmDouble>(pscalar); 375 } 306 376 307 377 /*create DoubleParam: */ … … 455 525 456 526 extern int my_rank; 527 extern int num_procs; 528 457 529 458 530 /*output: */ … … 473 545 #endif 474 546 547 /*cast to bool: */ 475 548 /*Assign output pointers: */ 476 549 *pboolean=(bool)booleanint; … … 482 555 483 556 extern int my_rank; 557 extern int num_procs; 484 558 485 559 /*output: */ … … 508 582 void IoModel::FetchData(IssmDouble* pscalar,int data_enum){ 509 583 584 510 585 extern int my_rank; 586 extern int num_procs; 587 511 588 512 589 /*output: */ 513 IssmPDouble scalar;514 int 590 IssmPDouble scalar; 591 int code; 515 592 516 593 /*Set file pointer to beginning of the data: */ … … 536 613 537 614 extern int my_rank; 615 extern int num_procs; 616 538 617 539 618 /*output: */ … … 547 626 if(code!=4)_error_("expecting a string for enum " << EnumToStringx(data_enum)); 548 627 628 /*Now fetch: */ 629 549 630 /*We have to read a string from disk. First read the dimensions of the string, then the string: */ 550 631 if(my_rank==0){ … … 574 655 } 575 656 657 576 658 /*Assign output pointers: */ 577 659 *pstring=string; … … 582 664 583 665 extern int my_rank; 666 extern int num_procs; 584 667 int i,j; 585 668 … … 597 680 if((code!=5) && (code!=6) && (code!=7))_error_("expecting a IssmDouble, integer or boolean matrix for enum " << EnumToStringx(data_enum)); 598 681 682 /*Now fetch: */ 683 599 684 /*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */ 600 685 /*numberofelements: */ … … 654 739 655 740 extern int my_rank; 741 extern int num_procs; 656 742 657 743 /*output: */ … … 665 751 if((code!=5) && (code!=6) && (code!=7))_error_("expecting a IssmDouble, integer or boolean matrix for enum " << EnumToStringx(data_enum)); 666 752 753 /*Now fetch: */ 754 667 755 /*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */ 668 756 /*numberofelements: */ … … 716 804 717 805 extern int my_rank; 806 extern int num_procs; 807 718 808 int i; 719 809 … … 784 874 785 875 int i; 876 786 877 extern int my_rank; 878 extern int num_procs; 787 879 788 880 /*output: */ 789 IssmDouble **matrices =NULL;790 int *mdims =NULL;791 int *ndims =NULL;792 int numrecords =0;881 IssmDouble** matrices=NULL; 882 int* mdims=NULL; 883 int* ndims=NULL; 884 int numrecords=0; 793 885 794 886 /*intermediary: */ … … 873 965 void IoModel::FetchData(Option** poption,int index){ 874 966 967 extern int my_rank; 968 extern int num_procs; 969 875 970 /*output: */ 876 int code;877 char *name= NULL;971 int code; 972 char *name = NULL; 878 973 879 974 /*First get option name*/ … … 919 1014 void IoModel::FetchData(int num,...){ 920 1015 921 va_list 922 int 923 IssmDouble *matrix =NULL;924 int 925 926 /*Go through the entire list of enums and fetch the corresponding data. 927 *Add it to the iomodel->data dataset. Everything928 * 1016 va_list ap; 1017 int dataenum; 1018 IssmDouble* matrix=NULL; 1019 int M,N; 1020 int i; 1021 1022 /*Go through the entire list of enums and fetch the corresponding data. Add it to the iomodel->data dataset. Everything 1023 *we fetch is a IssmDouble* : */ 929 1024 930 1025 va_start(ap,num); 931 for(i nt i=0; i<num; i++){1026 for(i=0; i<num; i++){ 932 1027 933 1028 dataenum=va_arg(ap, int); … … 953 1048 } 954 1049 va_end(ap); 1050 955 1051 } 956 1052 /*}}}*/ … … 970 1066 int nel; 971 1067 int numberofelements; 1068 972 1069 973 1070 /*variables being fetched: */ … … 1190 1287 1191 1288 extern int my_rank; 1289 extern int num_procs; 1192 1290 1193 1291 int found=0; … … 1250 1348 } 1251 1349 /*}}}*/ 1252 /*FUNCTION IoModel::DeclareIndependents{{{*/1253 void IoModel::DeclareIndependents(void){1254 1255 bool autodiff=false;1256 int dummy;1257 int i;1258 int num_independents;1259 int* independents=NULL;1260 1261 #ifdef _HAVE_ADOLC_1262 /*recover independent enums: */1263 this->Constant(&num_independents,AutodiffNumIndependentsEnum);1264 if(num_independents){1265 this->FetchData(&independents,&dummy,&dummy,AutodiffIndependentsEnum);1266 1267 /*now go fetch the independent variables for each independent enum: */1268 for(i=0;i<num_independents;i++){1269 this->FetchIndependent(independents[i]);1270 }1271 xDelete<int>(independents);1272 }1273 #else1274 /*if we asked for AD computations, we have a problem!: */1275 this->Constant(&autodiff,AutodiffIsautodiffEnum);1276 if(autodiff)_error_("Cannot carry out AD mode computations without support of ADOLC compiled in!");1277 #endif1278 1279 }1280 /*}}}*/1281 /*FUNCTION IoModel::FetchIndependent{{{*/1282 void IoModel::FetchIndependent(int independent_enum){1283 1284 #ifdef _HAVE_ADOLC_ //cannot come here unless you are running AD mode, from DeclaredIndependents:1285 extern int my_rank;1286 1287 /*output: */1288 int M,N;1289 IssmPDouble* buffer=NULL; //a buffer to read the data from disk1290 IssmDouble* matrix=NULL; //our independent variable1291 int code=0;1292 int vector_type=0;1293 1294 /*Set file pointer to beginning of the data: */1295 fid=this->SetFilePointerToData(&code,&vector_type,independent_enum);1296 if((code!=5) && (code!=6) && (code!=7))_error_("expecting a IssmDouble, integer or boolean matrix for enum " << EnumToStringx(independent_enum));1297 1298 /*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */1299 /*numberofelements: */1300 if(my_rank==0){1301 if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows for matrix ");1302 }1303 #ifdef _HAVE_MPI_1304 MPI_Bcast(&M,1,MPI_INT,0,MPI_COMM_WORLD);1305 #endif1306 1307 if(my_rank==0){1308 if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns for matrix ");1309 }1310 #ifdef _HAVE_MPI_1311 MPI_Bcast(&N,1,MPI_INT,0,MPI_COMM_WORLD);1312 #endif1313 1314 /*Now allocate matrix: */1315 if(M*N){1316 buffer=xNew<IssmPDouble>(M*N);1317 matrix=xNew<IssmDouble>(M*N);1318 1319 /*Read matrix on node 0, then broadcast: */1320 if(my_rank==0){1321 if(fread(buffer,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix ");1322 1323 /*Now, before we even broadcast this to other nodes, declare the whole matrix as a independent variable!: */1324 for (int i=0;i<M*N;++i) matrix[i]<<=buffer[i]; /*we use the <<= ADOLC overloaded operator to declare the independency*/1325 }1326 #ifdef _HAVE_MPI_1327 MPI_Bcast(matrix,M*N,MPI_DOUBLE,0,MPI_COMM_WORLD);1328 #endif1329 1330 xDelete<IssmPDouble>(buffer);1331 }1332 else _error_("cannot declare the independent variable " << EnumToStringx(independent_enum) << "if it's empty!");1333 1334 /*Ok, we are almost done. Matrix is now a independent matrix. We don't want this matrix to be fetched again in the1335 *future, which would effectively write over the independency in the ADOLC tape! So we are going to keep track of this1336 independent matrix inthe iomodel->data[independent_enum] data slot: */1337 this->data[independent_enum]=matrix;1338 this->independents[independent_enum]=true;1339 #endif1340 }1341 /*}}}*/ -
issm/trunk-jpl/src/c/classes/IoModel.h
r13413 r13432 39 39 /*for AD mode: to keep track of our independent variables we fetch:*/ 40 40 bool* independents; 41 DataSet* independent_objects; 41 42 42 43 /*Methods*/ -
issm/trunk-jpl/src/c/modules/AutodiffDriversx/AutodiffDriversx.cpp
r13362 r13432 12 12 13 13 void AutodiffDriversx(Elements* elements,Nodes* nodes,Vertices* vertices,Loads* loads,Materials* materials,Parameters* parameters,Results* results){ 14 bool isautodiff = false; 14 15 16 int i; 17 int dummy; 18 bool isautodiff = false; 19 int num_dependents; 20 int num_independents; 21 IssmDouble *axp = NULL; 22 int configuration_type; 23 int solveSize; 24 int edf_n ,edf_m; 25 double *xp = NULL; 26 int anIndepNum; 27 28 15 29 /*AD mode on?: */ 16 30 parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum); 31 17 32 if(isautodiff){ 33 18 34 #ifdef _HAVE_ADOLC_ 19 int num_dependents; 35 20 36 parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum); 21 int num_independents;22 37 parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum); 38 23 39 if(!(num_dependents*num_independents)) return; 24 int numberofvertices; 25 parameters->FindParam(&numberofvertices,MeshNumberofverticesEnum); 40 26 41 /*retrieve state variable: */ 27 IssmDouble *axp = NULL;28 parameters->FindParam(&axp,&num_independents,AutodiffXpEnum); 42 parameters->FindParam(&axp,&dummy,AutodiffXpEnum); 43 29 44 /* driver argument */ 30 double *xp=xNew<double>(num_independents);31 for(i nt i=0;i<num_independents;i++){45 xp=xNew<double>(num_independents); 46 for(i=0;i<num_independents;i++){ 32 47 xp[i]=reCast<double,IssmDouble>(axp[i]); 33 48 } 49 34 50 /* get the dimension for the solverx arguments*/ 35 int configuration_type;36 51 parameters->FindParam(&configuration_type,ConfigurationTypeEnum); 37 int solveSize=nodes->NumberOfDofs(configuration_type,FsetEnum); 38 int edf_n=solveSize*(solveSize+1), edf_m=solveSize; 52 solveSize=nodes->NumberOfDofs(configuration_type,FsetEnum); 53 edf_n=solveSize*(solveSize+1); 54 edf_m=solveSize; 39 55 40 56 /*get the EDF pointer:*/ … … 47 63 anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx; 48 64 // anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx; 65 49 66 /*allocate the space for the parameters to invoke the forward methods:*/ 50 67 anEDF_for_solverx_p->dp_x=xNew<double>(edf_n); … … 58 75 anEDF_for_solverx_p->dp_Z=xNew<double>(edf_n); 59 76 anEDF_for_solverx_p->dpp_Z=xNew<double>(num_dependents,edf_n); 77 60 78 /* Call AD driver:*/ 61 79 // single direction: 62 80 double *tangentDir=xNewZeroInit<double>(num_independents); 63 unsigned int anIndepNum=55; // <---------- make this selectable via config81 parameters->FindParam(&anIndepNum,AutodiffFosForwardIndexEnum); 64 82 tangentDir[anIndepNum]=1.0; 65 83 double *theJacVecProduct=xNew<double>(num_dependents); … … 67 85 if (fos_forward(1,num_dependents,num_independents, 0, xp, tangentDir, theOutput, theJacVecProduct )) 68 86 _error_("fos_forward returned non-zero error code"); 87 69 88 results->AddObject(new GenericExternalResult<IssmPDouble*>(results->Size()+1,AutodiffJacobianEnum,theJacVecProduct,num_dependents,1,1,0.0)); 70 89 … … 87 106 xDelete(anEDF_for_solverx_p->dp_Z); 88 107 xDelete(anEDF_for_solverx_p->dpp_Z); 108 89 109 /*Free ressources: */ 90 110 xDelete(xp); 91 xDelete(axp); // did we allocate this? 111 xDelete(axp); 112 92 113 #else 93 114 _error_("Should not be requesting AD drivers when an AD library is not available!"); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp
r13283 r13432 14 14 void CreateParametersAutodiff(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type){ 15 15 16 int i ,j;16 int i; 17 17 Parameters *parameters = NULL; 18 18 bool autodiff_analysis; 19 int* dependents = NULL; 20 int num_dependents; 21 int* independents = NULL; 22 int num_independents; 23 int numberofvertices; 19 int num_dependent_objects; 20 int num_dep=0; 21 int* names=NULL; 22 int* types=NULL; 23 int dummy; 24 24 25 IssmDouble* xp=NULL; 26 IssmDouble* xp_backup=NULL; 27 int num_ind,local_num_ind; 28 DataSet* dependent_objects=NULL; 25 29 26 30 /*Get parameters: */ … … 32 36 if(autodiff_analysis){ 33 37 34 iomodel->Constant(&num_independents,AutodiffNumIndependentsEnum);35 iomodel->Constant(&num_dependents,AutodiffNumDependentsEnum);36 iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum);38 /*retrieve driver: */ 39 parameters->AddObject(iomodel->CopyConstantObject(AutodiffDriverEnum)); 40 parameters->AddObject(iomodel->CopyConstantObject(AutodiffFosForwardIndexEnum)); 37 41 38 /*recover dependents: */ 39 parameters->AddObject(iomodel->CopyConstantObject(AutodiffNumDependentsEnum)); 40 if(num_dependents){ 41 iomodel->FetchData(&dependents,NULL,&num_dependents,AutodiffDependentsEnum); 42 parameters->AddObject(new IntVecParam(AutodiffDependentsEnum,dependents,num_dependents)); 42 /*Deal with dependents first: {{{*/ 43 iomodel->Constant(&num_dependent_objects,AutodiffNumDependentObjectsEnum); 44 dependent_objects=new DataSet(); 45 num_dep=0; 46 47 if(num_dependent_objects){ 48 iomodel->FetchData(&names,&dummy,&dummy,AutodiffDependentObjectNamesEnum); 49 iomodel->FetchData(&types,&dummy,&dummy,AutodiffDependentObjectTypesEnum); 50 51 for(i=0;i<num_dependent_objects;i++){ 52 DependentObject* dep=new DependentObject(names[i],types[i]); 53 dependent_objects->AddObject(dep); 54 num_dep+=dep->NumDependents(); 55 } 56 57 /*Free ressources:*/ 58 xDelete<int>(names); 59 xDelete<int>(types); 43 60 } 61 parameters->AddObject(new DataSetParam(AutodiffDependentObjectsEnum,dependent_objects)); 62 parameters->AddObject(new IntParam(AutodiffNumDependentsEnum,num_dep)); 44 63 45 /*recover independents: */ 46 parameters->AddObject(iomodel->CopyConstantObject(AutodiffNumIndependentsEnum)); 47 if(num_independents){ 48 iomodel->FetchData(&independents,NULL,&num_independents,AutodiffIndependentsEnum); 49 parameters->AddObject(new IntVecParam(AutodiffIndependentsEnum,independents,num_independents)); 64 delete dependent_objects; 65 /*}}}*/ 50 66 51 /*Build state vector, value at which we compute our gradients of dependent variables in adolc: the xp vector */ 52 xp=xNew<IssmDouble>(num_independents*numberofvertices); 53 for(i=0;i<num_independents;i++){ 54 IssmDouble* values=iomodel->data[independents[i]]; 55 for(j=0;j<numberofvertices;j++){ 56 xp[i*numberofvertices+j]=values[j]; 57 } 67 /*Deal with independents: {{{*/ 68 69 /*Independents have already been recovered in iomodel->DeclareIndependents. Just do some more processing. 70 *In particular, figure out num_independents, and create the state vector xp, or size num_independents x 1 :*/ 71 num_ind=0; 72 for(i=0;i<iomodel->independent_objects->Size();i++){ 73 IndependentObject* ind=(IndependentObject*)iomodel->independent_objects->GetObjectByOffset(i); 74 num_ind+=ind->NumIndependents(); 75 } 76 if(num_ind){ 77 xp=xNew<IssmDouble>(num_ind); 78 xp_backup=xp; 79 for(i=0;i<iomodel->independent_objects->Size();i++){ 80 IndependentObject* ind=(IndependentObject*)iomodel->independent_objects->GetObjectByOffset(i); 81 ind->FillIndependents(iomodel->data,xp); 82 local_num_ind=ind->NumIndependents(); xp=xp+local_num_ind; 58 83 } 59 parameters->AddObject(new DoubleVecParam(AutodiffXpEnum,xp,num_independents*numberofvertices));84 xp=xp_backup; parameters->AddObject(new DoubleVecParam(AutodiffXpEnum,xp,num_ind)); 60 85 } 86 parameters->AddObject(new IntParam(AutodiffNumIndependentsEnum,num_ind)); 87 88 /*Don't forget to copy iomodel->independent_objects to parameters: */ 89 parameters->AddObject(new DataSetParam(AutodiffIndependentObjectsEnum,iomodel->independent_objects)); 90 /*}}}*/ 61 91 62 92 /*Assign output pointer: */ -
issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp
r13283 r13432 20 20 21 21 int i,analysis_type,dim,verbose; 22 bool isthermal,isprognostic,isdiagnostic,isgroundingline,isenthalpy ,autodiff;22 bool isthermal,isprognostic,isdiagnostic,isgroundingline,isenthalpy; 23 23 24 24 /*output: */ … … 43 43 iomodel->Constant(&isdiagnostic,TransientIsdiagnosticEnum); 44 44 iomodel->Constant(&isgroundingline,TransientIsgroundinglineEnum); 45 iomodel->Constant(&autodiff,AutodiffIsautodiffEnum);46 47 /*If we are running in AD mode, we need to declare our independent variables now, before48 *and prevent them from being erased during successive calls to iomodel->FetchData and49 iomodel->DeleteData:*/50 if(autodiff)iomodel->DeclareIndependents();51 45 52 46 SetVerbosityLevel(verbose); -
issm/trunk-jpl/src/c/modules/RequestedDependentsx/RequestedDependentsx.cpp
r13285 r13432 14 14 15 15 int i; 16 int output_enum;17 16 bool isautodiff = false; 18 17 IssmDouble output_value; 19 Element *element = NULL;20 int *dependent_enums = NULL;21 18 22 19 int num_dependents; 23 20 IssmPDouble *dependents; 21 DataSet* dependent_objects=NULL; 24 22 25 23 /*AD mode on?: */ … … 29 27 #ifdef _HAVE_ADOLC_ 30 28 parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum); 29 parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum); 31 30 if(num_dependents){ 32 31 dependents=xNew<IssmPDouble>(num_dependents); 33 parameters->FindParam(&dependent_enums,&num_dependents,AutodiffDependentsEnum);34 32 35 33 /*Go through our dependent variables, and compute the response:*/ 36 for(i=0;i<num_dependents;i++){ 37 Responsex(&output_value,elements,nodes,vertices,loads,materials,parameters,dependent_enums[i],false,0); 34 for(i=0;i<dependent_objects->Size();i++){ 35 DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i); 36 Responsex(&output_value,elements,nodes,vertices,loads,materials,parameters,dep->name,false,0); 38 37 output_value>>=dependents[i]; 39 38 } 40 39 } 40 delete dependent_objects; 41 41 #else 42 42 _error_("Should not be requesting dependents when an AD library is not available!");
Note:
See TracChangeset
for help on using the changeset viewer.
