Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 13431) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 13432) @@ -54,4 +54,8 @@ ./classes/objects/Constraints/SpcTransient.cpp\ ./classes/objects/Constraints/SpcTransient.h\ + ./classes/objects/IndependentObject.h\ + ./classes/objects/IndependentObject.cpp\ + ./classes/objects/DependentObject.h\ + ./classes/objects/DependentObject.cpp\ ./classes/DofIndexing.h\ ./classes/DofIndexing.cpp\ @@ -152,4 +156,6 @@ ./classes/objects/Params/TransientParam.h\ ./classes/objects/Params/TransientParam.cpp\ + ./classes/objects/Params/DataSetParam.h\ + ./classes/objects/Params/DataSetParam.cpp\ ./Container/Container.h\ ./Container/Constraints.h\ Index: /issm/trunk-jpl/src/c/classes/IoModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 13431) +++ /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 13432) @@ -24,23 +24,23 @@ /*FUNCTION IoModel::IoModel(){{{*/ IoModel::IoModel(){ - this->fid = NULL; - this->data = NULL; - this->independents = NULL; - this->constants = NULL; - - this->my_elements = NULL; - this->my_nodes = NULL; - this->my_vertices = NULL; - this->singlenodetoelementconnectivity = NULL; - this->numbernodetoelementconnectivity = NULL; - - this->nodecounter = 0; - this->loadcounter = 0; - this->constraintcounter = 0; + this->fid=NULL; + this->data=NULL; + this->independents=NULL; + this->constants=NULL; + + this->my_elements=NULL; + this->my_nodes=NULL; + this->my_vertices=NULL; + this->singlenodetoelementconnectivity=NULL; + this->numbernodetoelementconnectivity=NULL; + + this->nodecounter=0; + this->loadcounter=0; + this->constraintcounter=0; } /*}}}*/ /*FUNCTION IoModel::IoModel(FILE* iomodel_handle){{{*/ IoModel::IoModel(FILE* iomodel_handle){ - + /*First, keep track of the file handle: */ this->fid=iomodel_handle; @@ -48,4 +48,13 @@ /*Check that Enums are Synchronized*/ this->CheckEnumSync(); + + /*Initialize data: */ + this->data=xNew(MaximumNumberOfEnums); + for(int i=0;idata[i]=NULL; + + /*If we are running in AD mode, we need to declare our independent variables now, + *and prevent them from being erased during successive calls to iomodel->FetchConstants, iomodel->FetchData and + iomodel->DeleteData:*/ + this->DeclareIndependents(); /*Initialize and read constants:*/ @@ -53,22 +62,14 @@ this->FetchConstants(); /*this routine goes through the input file, and fetches bools, ints, IssmDoubles and strings only, nothing memory intensive*/ - /*Initialize data: */ - this->data=xNew(MaximumNumberOfEnums); - for(int i=0;idata[i]=NULL; - - /*Initialize array detecting whether data[i] is an independent AD mode variable: */ - this->independents=xNew(MaximumNumberOfEnums); - for(int i=0;iindependents[i]=false; - /*Initialize permanent data: */ - this->my_elements = NULL; - this->my_nodes = NULL; - this->my_vertices = NULL; - this->singlenodetoelementconnectivity = NULL; - this->numbernodetoelementconnectivity = NULL; - - this->nodecounter = 0; - this->loadcounter = 0; - this->constraintcounter = 0; + this->my_elements=NULL; + this->my_nodes=NULL; + this->my_vertices=NULL; + this->singlenodetoelementconnectivity=NULL; + this->numbernodetoelementconnectivity=NULL; + + this->nodecounter=0; + this->loadcounter=0; + this->constraintcounter=0; } /*}}}*/ @@ -91,4 +92,5 @@ xDelete(this->data); xDelete(this->independents); + delete this->independent_objects; xDelete(this->my_elements); xDelete(this->my_nodes); @@ -197,14 +199,72 @@ } /*}}}*/ +/*FUNCTION IoModel::DeclareIndependents{{{*/ +void IoModel::DeclareIndependents(void){ + + + int i; + bool autodiff = false; + int num_independent_objects; + + int *names = NULL; + int *types = NULL; + + int numberofvertices; + int dummy; + + + /*Initialize array detecting whether data[i] is an independent AD mode variable: */ + this->independents=xNew(MaximumNumberOfEnums); + for(i=0;iindependents[i]=false; + + this->FetchData(&autodiff,AutodiffIsautodiffEnum); + if(autodiff){ + + this->FetchData(&numberofvertices,MeshNumberofverticesEnum); + + #ifdef _HAVE_ADOLC_ + /*build dataset made of independent objects: */ + this->independent_objects=new DataSet(); + this->FetchData(&num_independent_objects,AutodiffNumIndependentObjectsEnum); + if(num_independent_objects){ + this->FetchData(&names,&dummy,&dummy,AutodiffIndependentObjectNamesEnum); + this->FetchData(&types,&dummy,&dummy,AutodiffIndependentObjectTypesEnum); + + /*create independent objects, and at the same time, fetch the corresponding independent variables, + *and declare them as such in ADOLC: */ + for(i=0;iindependent_objects->AddObject(independent_object); + + /*now go fetch the independent variable: */ + independent_object->FetchIndependent(this); //supply the pointer to iomodel. + } + xDelete(names); + xDelete(types); + } + #else + /*if we asked for AD computations, we have a problem!: */ + _error_("Cannot carry out AD mode computations without support of ADOLC compiled in!"); + #endif + } + +} +/*}}}*/ /*FUNCTION IoModel::DeleteData(int num,...){{{*/ void IoModel::DeleteData(int num,...){ - va_list ap; - int dataenum; - DoubleMatParam *parameter = NULL; + va_list ap; + int dataenum; + int i; + DoubleMatParam* parameter=NULL; /*Go through the entire list of enums and delete the corresponding data from the iomodel-data dataset: */ + va_start(ap,num); - for(int i=0; ifid)!=1) _error_("could not read scalar "); - #ifdef _HAVE_MPI_ - MPI_Bcast(&scalar,1,MPI_DOUBLE,0,MPI_COMM_WORLD); - #endif + /*Read the scalar and broadcast it to other cpus. However, if this record has already + * been read in DeclareIndependents, we should not re-read it, but grab it from the iomodel->data + * slots: */ + if(this->independents[record_enum]){ + scalar=*(this->data[record_enum]); + } + else{ + if(fread(&pscalar,sizeof(IssmPDouble),1,this->fid)!=1) _error_("could not read scalar "); + #ifdef _HAVE_MPI_ + MPI_Bcast(&pscalar,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + #endif + scalar=reCast(pscalar); + } /*create DoubleParam: */ @@ -455,4 +525,6 @@ extern int my_rank; + extern int num_procs; + /*output: */ @@ -473,4 +545,5 @@ #endif + /*cast to bool: */ /*Assign output pointers: */ *pboolean=(bool)booleanint; @@ -482,4 +555,5 @@ extern int my_rank; + extern int num_procs; /*output: */ @@ -508,9 +582,12 @@ void IoModel::FetchData(IssmDouble* pscalar,int data_enum){ + extern int my_rank; + extern int num_procs; + /*output: */ - IssmPDouble scalar; - int code; + IssmPDouble scalar; + int code; /*Set file pointer to beginning of the data: */ @@ -536,4 +613,6 @@ extern int my_rank; + extern int num_procs; + /*output: */ @@ -547,4 +626,6 @@ if(code!=4)_error_("expecting a string for enum " << EnumToStringx(data_enum)); + /*Now fetch: */ + /*We have to read a string from disk. First read the dimensions of the string, then the string: */ if(my_rank==0){ @@ -574,4 +655,5 @@ } + /*Assign output pointers: */ *pstring=string; @@ -582,4 +664,5 @@ extern int my_rank; + extern int num_procs; int i,j; @@ -597,4 +680,6 @@ if((code!=5) && (code!=6) && (code!=7))_error_("expecting a IssmDouble, integer or boolean matrix for enum " << EnumToStringx(data_enum)); + /*Now fetch: */ + /*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */ /*numberofelements: */ @@ -654,4 +739,5 @@ extern int my_rank; + extern int num_procs; /*output: */ @@ -665,4 +751,6 @@ if((code!=5) && (code!=6) && (code!=7))_error_("expecting a IssmDouble, integer or boolean matrix for enum " << EnumToStringx(data_enum)); + /*Now fetch: */ + /*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */ /*numberofelements: */ @@ -716,4 +804,6 @@ extern int my_rank; + extern int num_procs; + int i; @@ -784,11 +874,13 @@ int i; + extern int my_rank; + extern int num_procs; /*output: */ - IssmDouble **matrices = NULL; - int *mdims = NULL; - int *ndims = NULL; - int numrecords = 0; + IssmDouble** matrices=NULL; + int* mdims=NULL; + int* ndims=NULL; + int numrecords=0; /*intermediary: */ @@ -873,7 +965,10 @@ void IoModel::FetchData(Option** poption,int index){ + extern int my_rank; + extern int num_procs; + /*output: */ - int code; - char *name = NULL; + int code; + char *name = NULL; /*First get option name*/ @@ -919,15 +1014,15 @@ void IoModel::FetchData(int num,...){ - va_list ap; - int dataenum; - IssmDouble *matrix = NULL; - int M,N; - - /*Go through the entire list of enums and fetch the corresponding data. - * Add it to the iomodel->data dataset. Everything - * we fetch is a IssmDouble* : */ + va_list ap; + int dataenum; + IssmDouble* matrix=NULL; + int M,N; + int i; + + /*Go through the entire list of enums and fetch the corresponding data. Add it to the iomodel->data dataset. Everything + *we fetch is a IssmDouble* : */ va_start(ap,num); - for(int i=0; iConstant(&num_independents,AutodiffNumIndependentsEnum); - if(num_independents){ - this->FetchData(&independents,&dummy,&dummy,AutodiffIndependentsEnum); - - /*now go fetch the independent variables for each independent enum: */ - for(i=0;iFetchIndependent(independents[i]); - } - xDelete(independents); - } - #else - /*if we asked for AD computations, we have a problem!: */ - this->Constant(&autodiff,AutodiffIsautodiffEnum); - if(autodiff)_error_("Cannot carry out AD mode computations without support of ADOLC compiled in!"); - #endif - -} -/*}}}*/ -/*FUNCTION IoModel::FetchIndependent{{{*/ -void IoModel::FetchIndependent(int independent_enum){ - - #ifdef _HAVE_ADOLC_ //cannot come here unless you are running AD mode, from DeclaredIndependents: - extern int my_rank; - - /*output: */ - int M,N; - IssmPDouble* buffer=NULL; //a buffer to read the data from disk - IssmDouble* matrix=NULL; //our independent variable - int code=0; - int vector_type=0; - - /*Set file pointer to beginning of the data: */ - fid=this->SetFilePointerToData(&code,&vector_type,independent_enum); - if((code!=5) && (code!=6) && (code!=7))_error_("expecting a IssmDouble, integer or boolean matrix for enum " << EnumToStringx(independent_enum)); - - /*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */ - /*numberofelements: */ - if(my_rank==0){ - if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows for matrix "); - } - #ifdef _HAVE_MPI_ - MPI_Bcast(&M,1,MPI_INT,0,MPI_COMM_WORLD); - #endif - - if(my_rank==0){ - if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns for matrix "); - } - #ifdef _HAVE_MPI_ - MPI_Bcast(&N,1,MPI_INT,0,MPI_COMM_WORLD); - #endif - - /*Now allocate matrix: */ - if(M*N){ - buffer=xNew(M*N); - matrix=xNew(M*N); - - /*Read matrix on node 0, then broadcast: */ - if(my_rank==0){ - if(fread(buffer,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix "); - - /*Now, before we even broadcast this to other nodes, declare the whole matrix as a independent variable!: */ - for (int i=0;i(buffer); - } - else _error_("cannot declare the independent variable " << EnumToStringx(independent_enum) << "if it's empty!"); - - /*Ok, we are almost done. Matrix is now a independent matrix. We don't want this matrix to be fetched again in the - *future, which would effectively write over the independency in the ADOLC tape! So we are going to keep track of this - independent matrix inthe iomodel->data[independent_enum] data slot: */ - this->data[independent_enum]=matrix; - this->independents[independent_enum]=true; - #endif -} -/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/IoModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/IoModel.h (revision 13431) +++ /issm/trunk-jpl/src/c/classes/IoModel.h (revision 13432) @@ -39,4 +39,5 @@ /*for AD mode: to keep track of our independent variables we fetch:*/ bool* independents; + DataSet* independent_objects; /*Methods*/ Index: /issm/trunk-jpl/src/c/modules/AutodiffDriversx/AutodiffDriversx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/AutodiffDriversx/AutodiffDriversx.cpp (revision 13431) +++ /issm/trunk-jpl/src/c/modules/AutodiffDriversx/AutodiffDriversx.cpp (revision 13432) @@ -12,29 +12,45 @@ void AutodiffDriversx(Elements* elements,Nodes* nodes,Vertices* vertices,Loads* loads,Materials* materials,Parameters* parameters,Results* results){ - bool isautodiff = false; + + + int i; + int dummy; + bool isautodiff = false; + int num_dependents; + int num_independents; + IssmDouble *axp = NULL; + int configuration_type; + int solveSize; + int edf_n ,edf_m; + double *xp = NULL; + int anIndepNum; + + /*AD mode on?: */ parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum); + if(isautodiff){ + #ifdef _HAVE_ADOLC_ - int num_dependents; + parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum); - int num_independents; parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum); + if(!(num_dependents*num_independents)) return; - int numberofvertices; - parameters->FindParam(&numberofvertices,MeshNumberofverticesEnum); + /*retrieve state variable: */ - IssmDouble *axp = NULL; - parameters->FindParam(&axp,&num_independents,AutodiffXpEnum); + parameters->FindParam(&axp,&dummy,AutodiffXpEnum); + /* driver argument */ - double *xp=xNew(num_independents); - for(int i=0;i(num_independents); + for(i=0;i(axp[i]); } + /* get the dimension for the solverx arguments*/ - int configuration_type; parameters->FindParam(&configuration_type,ConfigurationTypeEnum); - int solveSize=nodes->NumberOfDofs(configuration_type,FsetEnum); - int edf_n=solveSize*(solveSize+1), edf_m=solveSize; + solveSize=nodes->NumberOfDofs(configuration_type,FsetEnum); + edf_n=solveSize*(solveSize+1); + edf_m=solveSize; /*get the EDF pointer:*/ @@ -47,4 +63,5 @@ anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx; // anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx; + /*allocate the space for the parameters to invoke the forward methods:*/ anEDF_for_solverx_p->dp_x=xNew(edf_n); @@ -58,8 +75,9 @@ anEDF_for_solverx_p->dp_Z=xNew(edf_n); anEDF_for_solverx_p->dpp_Z=xNew(num_dependents,edf_n); + /* Call AD driver:*/ // single direction: double *tangentDir=xNewZeroInit(num_independents); - unsigned int anIndepNum=55; // <---------- make this selectable via config + parameters->FindParam(&anIndepNum,AutodiffFosForwardIndexEnum); tangentDir[anIndepNum]=1.0; double *theJacVecProduct=xNew(num_dependents); @@ -67,4 +85,5 @@ if (fos_forward(1,num_dependents,num_independents, 0, xp, tangentDir, theOutput, theJacVecProduct )) _error_("fos_forward returned non-zero error code"); + results->AddObject(new GenericExternalResult(results->Size()+1,AutodiffJacobianEnum,theJacVecProduct,num_dependents,1,1,0.0)); @@ -87,7 +106,9 @@ xDelete(anEDF_for_solverx_p->dp_Z); xDelete(anEDF_for_solverx_p->dpp_Z); + /*Free ressources: */ xDelete(xp); - xDelete(axp); // did we allocate this? + xDelete(axp); + #else _error_("Should not be requesting AD drivers when an AD library is not available!"); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp (revision 13431) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp (revision 13432) @@ -14,13 +14,17 @@ void CreateParametersAutodiff(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type){ - int i,j; + int i; Parameters *parameters = NULL; bool autodiff_analysis; - int* dependents = NULL; - int num_dependents; - int* independents = NULL; - int num_independents; - int numberofvertices; + int num_dependent_objects; + int num_dep=0; + int* names=NULL; + int* types=NULL; + int dummy; + IssmDouble* xp=NULL; + IssmDouble* xp_backup=NULL; + int num_ind,local_num_ind; + DataSet* dependent_objects=NULL; /*Get parameters: */ @@ -32,31 +36,57 @@ if(autodiff_analysis){ - iomodel->Constant(&num_independents,AutodiffNumIndependentsEnum); - iomodel->Constant(&num_dependents,AutodiffNumDependentsEnum); - iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum); + /*retrieve driver: */ + parameters->AddObject(iomodel->CopyConstantObject(AutodiffDriverEnum)); + parameters->AddObject(iomodel->CopyConstantObject(AutodiffFosForwardIndexEnum)); - /*recover dependents: */ - parameters->AddObject(iomodel->CopyConstantObject(AutodiffNumDependentsEnum)); - if(num_dependents){ - iomodel->FetchData(&dependents,NULL,&num_dependents,AutodiffDependentsEnum); - parameters->AddObject(new IntVecParam(AutodiffDependentsEnum,dependents,num_dependents)); + /*Deal with dependents first: {{{*/ + iomodel->Constant(&num_dependent_objects,AutodiffNumDependentObjectsEnum); + dependent_objects=new DataSet(); + num_dep=0; + + if(num_dependent_objects){ + iomodel->FetchData(&names,&dummy,&dummy,AutodiffDependentObjectNamesEnum); + iomodel->FetchData(&types,&dummy,&dummy,AutodiffDependentObjectTypesEnum); + + for(i=0;iAddObject(dep); + num_dep+=dep->NumDependents(); + } + + /*Free ressources:*/ + xDelete(names); + xDelete(types); } + parameters->AddObject(new DataSetParam(AutodiffDependentObjectsEnum,dependent_objects)); + parameters->AddObject(new IntParam(AutodiffNumDependentsEnum,num_dep)); - /*recover independents: */ - parameters->AddObject(iomodel->CopyConstantObject(AutodiffNumIndependentsEnum)); - if(num_independents){ - iomodel->FetchData(&independents,NULL,&num_independents,AutodiffIndependentsEnum); - parameters->AddObject(new IntVecParam(AutodiffIndependentsEnum,independents,num_independents)); + delete dependent_objects; + /*}}}*/ - /*Build state vector, value at which we compute our gradients of dependent variables in adolc: the xp vector */ - xp=xNew(num_independents*numberofvertices); - for(i=0;idata[independents[i]]; - for(j=0;jDeclareIndependents. Just do some more processing. + *In particular, figure out num_independents, and create the state vector xp, or size num_independents x 1 :*/ + num_ind=0; + for(i=0;iindependent_objects->Size();i++){ + IndependentObject* ind=(IndependentObject*)iomodel->independent_objects->GetObjectByOffset(i); + num_ind+=ind->NumIndependents(); + } + if(num_ind){ + xp=xNew(num_ind); + xp_backup=xp; + for(i=0;iindependent_objects->Size();i++){ + IndependentObject* ind=(IndependentObject*)iomodel->independent_objects->GetObjectByOffset(i); + ind->FillIndependents(iomodel->data,xp); + local_num_ind=ind->NumIndependents(); xp=xp+local_num_ind; } - parameters->AddObject(new DoubleVecParam(AutodiffXpEnum,xp,num_independents*numberofvertices)); + xp=xp_backup; parameters->AddObject(new DoubleVecParam(AutodiffXpEnum,xp,num_ind)); } + parameters->AddObject(new IntParam(AutodiffNumIndependentsEnum,num_ind)); + + /*Don't forget to copy iomodel->independent_objects to parameters: */ + parameters->AddObject(new DataSetParam(AutodiffIndependentObjectsEnum,iomodel->independent_objects)); + /*}}}*/ /*Assign output pointer: */ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 13431) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 13432) @@ -20,5 +20,5 @@ int i,analysis_type,dim,verbose; - bool isthermal,isprognostic,isdiagnostic,isgroundingline,isenthalpy,autodiff; + bool isthermal,isprognostic,isdiagnostic,isgroundingline,isenthalpy; /*output: */ @@ -43,10 +43,4 @@ iomodel->Constant(&isdiagnostic,TransientIsdiagnosticEnum); iomodel->Constant(&isgroundingline,TransientIsgroundinglineEnum); - iomodel->Constant(&autodiff,AutodiffIsautodiffEnum); - - /*If we are running in AD mode, we need to declare our independent variables now, before - *and prevent them from being erased during successive calls to iomodel->FetchData and - iomodel->DeleteData:*/ - if(autodiff)iomodel->DeclareIndependents(); SetVerbosityLevel(verbose); Index: /issm/trunk-jpl/src/c/modules/RequestedDependentsx/RequestedDependentsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/RequestedDependentsx/RequestedDependentsx.cpp (revision 13431) +++ /issm/trunk-jpl/src/c/modules/RequestedDependentsx/RequestedDependentsx.cpp (revision 13432) @@ -14,12 +14,10 @@ int i; - int output_enum; bool isautodiff = false; IssmDouble output_value; - Element *element = NULL; - int *dependent_enums = NULL; int num_dependents; IssmPDouble *dependents; + DataSet* dependent_objects=NULL; /*AD mode on?: */ @@ -29,14 +27,16 @@ #ifdef _HAVE_ADOLC_ parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum); + parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum); if(num_dependents){ dependents=xNew(num_dependents); - parameters->FindParam(&dependent_enums,&num_dependents,AutodiffDependentsEnum); /*Go through our dependent variables, and compute the response:*/ - for(i=0;iSize();i++){ + DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i); + Responsex(&output_value,elements,nodes,vertices,loads,materials,parameters,dep->name,false,0); output_value>>=dependents[i]; } } + delete dependent_objects; #else _error_("Should not be requesting dependents when an AD library is not available!");