Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 18871) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 18872) @@ -43,9 +43,7 @@ /*Object constructors and destructor*/ -FemModel::FemModel(int argc,char** argv,ISSM_MPI_Comm incomm){/*{{{*/ +FemModel::FemModel(int argc,char** argv,ISSM_MPI_Comm incomm,bool trace){/*{{{*/ /*configuration: */ - int* analyses=NULL; - int numanalyses; int solution_type; int ierr; @@ -77,9 +75,8 @@ /*Create femmodel from input files: */ profiler->Tag(StartInit); - this->InitFromFiles(rootpath,binfilename,outbinfilename,petscfilename,lockfilename,solution_type); + this->InitFromFiles(rootpath,binfilename,outbinfilename,petscfilename,lockfilename,solution_type,trace,NULL); profiler->Tag(FinishInit); /*Free resources */ - xDelete(analyses); xDelete(lockfilename); xDelete(binfilename); @@ -87,4 +84,18 @@ xDelete(petscfilename); xDelete(rootpath); + +} +/*}}}*/ +FemModel::FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, ISSM_MPI_Comm incomm, int solution_type,IssmPDouble* X){ /*{{{*/ + + bool traceon=true; + + /*Store the communicator, but do not set it as a global variable, as this has already + * been done by the FemModel that called this copy constructor: */ + this->comm=incomm; + this->SetStaticComm(); + + /*Create femmodel from input files, with trace activated: */ + this->InitFromFiles(rootpath,inputfilename,outputfilename,toolkitsfilename,lockfilename,solution_type,traceon,X); } @@ -131,5 +142,5 @@ } /*}}}*/ -void FemModel::InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, const int in_solution_type){/*{{{*/ +void FemModel::InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, const int in_solution_type,bool trace,IssmPDouble* X){/*{{{*/ /*intermediary*/ @@ -148,8 +159,7 @@ this->analysis_counter = nummodels-1; //point to last analysis_type carried out. this->results = new Results(); //not initialized by CreateDataSets - /*Open input file on cpu 0 and create IoModel */ if(my_rank==0) IOMODEL = pfopen0(inputfilename ,"rb"); - IoModel* iomodel = new IoModel(IOMODEL); + IoModel* iomodel = new IoModel(IOMODEL,trace,X); /*Figure out what analyses are activated for this solution*/ @@ -188,10 +198,13 @@ pfclose(toolkitsoptionsfid,toolkitsfilename); - /*Open output file once for all and add output file name and file descriptor to parameters*/ + /*Open output file once for all and add output file descriptor to parameters*/ output_fid=pfopen(outputfilename,"wb"); + this->parameters->SetParam(output_fid,OutputFilePointerEnum); + + /*Now save all of these file names into parameters, you never know when you might need them: */ + this->parameters->AddObject(new StringParam(ToolkitsFileNameEnum,toolkitsfilename)); + this->parameters->AddObject(new StringParam(RootPathEnum,rootpath)); + this->parameters->AddObject(new StringParam(InputFileNameEnum,inputfilename)); this->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename)); - this->parameters->SetParam(output_fid,OutputFilePointerEnum); - - /*Save lock file name for later: */ this->parameters->AddObject(new StringParam(LockFileNameEnum,lockfilename)); @@ -741,24 +754,33 @@ } else{ - switch(output_enum){ - - /*Scalar output*/ - case DivergenceEnum: this->Divergencex(&double_result); break; - case MaxDivergenceEnum: this->MaxDivergencex(&double_result); break; - case IceVolumeEnum: this->IceVolumex(&double_result); break; - case IceVolumeAboveFloatationEnum: this->IceVolumeAboveFloatationx(&double_result); break; - case MinVelEnum: this->MinVelx(&double_result); break; - case MaxVelEnum: this->MaxVelx(&double_result); break; - case MinVxEnum: this->MinVxx(&double_result); break; - case MaxVxEnum: this->MaxVxx(&double_result); break; - case MaxAbsVxEnum: this->MaxAbsVxx(&double_result); break; - case MinVyEnum: this->MinVyx(&double_result); break; - case MaxVyEnum: this->MaxVyx(&double_result); break; - case MaxAbsVyEnum: this->MaxAbsVyx(&double_result); break; - case MinVzEnum: this->MinVzx(&double_result); break; - case MaxVzEnum: this->MaxVzx(&double_result); break; - case MaxAbsVzEnum: this->MaxAbsVzx(&double_result); break; - case MassFluxEnum: this->MassFluxx(&double_result); break; - case TotalSmbEnum: this->TotalSmbx(&double_result); break; + /*last chance for the output definition, if the enum is one of Outputdefinition[1-10]Enum:*/ + if(output_enum>=Outputdefinition1Enum && output_enum <=Outputdefinition10Enum){ + double_result = OutputDefinitionsResponsex(this,output_enum); + if(save_results){ + results->AddResult(new GenericExternalResult(results->Size()+1,output_string,reCast(double_result),step,time)); + continue; + } + } + else{ + switch(output_enum){ + + /*Scalar output*/ + case DivergenceEnum: this->Divergencex(&double_result); break; + case MaxDivergenceEnum: this->MaxDivergencex(&double_result); break; + case IceVolumeEnum: this->IceVolumex(&double_result); break; + case IceVolumeAboveFloatationEnum: this->IceVolumeAboveFloatationx(&double_result); break; + case MinVelEnum: this->MinVelx(&double_result); break; + case MaxVelEnum: this->MaxVelx(&double_result); break; + case MinVxEnum: this->MinVxx(&double_result); break; + case MaxVxEnum: this->MaxVxx(&double_result); break; + case MaxAbsVxEnum: this->MaxAbsVxx(&double_result); break; + case MinVyEnum: this->MinVyx(&double_result); break; + case MaxVyEnum: this->MaxVyx(&double_result); break; + case MaxAbsVyEnum: this->MaxAbsVyx(&double_result); break; + case MinVzEnum: this->MinVzx(&double_result); break; + case MaxVzEnum: this->MaxVzx(&double_result); break; + case MaxAbsVzEnum: this->MaxAbsVzx(&double_result); break; + case MassFluxEnum: this->MassFluxx(&double_result); break; + case TotalSmbEnum: this->TotalSmbx(&double_result); break; /*Scalar control output*/ @@ -778,67 +800,69 @@ case SurfaceAbsMisfitEnum: SurfaceAbsMisfitx(&double_result); break; - /*Vector */ - default: - - /*Vector layout*/ - int interpolation,nodesperelement,size; - int rank_interpolation=-1,rank_nodesperelement=-1; - - /*Get interpolation (and compute input if necessary)*/ - for(int j=0;jSize();j++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(j)); - element->ResultInterpolation(&rank_interpolation,&rank_nodesperelement,output_enum); - } - - /*Broadcast for cpus that do not have any elements*/ - ISSM_MPI_Reduce(&rank_interpolation,&interpolation,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm()); - ISSM_MPI_Reduce(&rank_nodesperelement,&nodesperelement,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(&interpolation,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(&nodesperelement,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - - if(results_on_nodes){ - - /*Allocate matrices*/ - int nbe = this->elements->NumberOfElements(); - IssmDouble* values = xNewZeroInit(nbe*nodesperelement); - IssmDouble* allvalues = xNew(nbe*nodesperelement); - - /*Fill-in matrix*/ + /*Vector */ + default: + + /*Vector layout*/ + int interpolation,nodesperelement,size; + int rank_interpolation=-1,rank_nodesperelement=-1; + + /*Get interpolation (and compute input if necessary)*/ for(int j=0;jSize();j++){ Element* element=xDynamicCast(this->elements->GetObjectByOffset(j)); - element->ResultToPatch(values,nodesperelement,output_enum); + element->ResultInterpolation(&rank_interpolation,&rank_nodesperelement,output_enum); } - /*Gather from all cpus*/ - ISSM_MPI_Allreduce((void*)values,(void*)allvalues,nbe*nodesperelement,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); - xDelete(values); - - if(save_results)results->AddResult(new GenericExternalResult(results->Size()+1,output_enum,allvalues,nbe,nodesperelement,step,time)); - xDelete(allvalues); - - } - else{ - - /*Allocate vector depending on interpolation*/ - switch(interpolation){ - case P0Enum: size = this->elements->NumberOfElements(); break; - case P1Enum: size = this->vertices->NumberOfVertices(); break; - default: _error_("Interpolation "<elements->NumberOfElements(); + IssmDouble* values = xNewZeroInit(nbe*nodesperelement); + IssmDouble* allvalues = xNew(nbe*nodesperelement); + + /*Fill-in matrix*/ + for(int j=0;jSize();j++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(j)); + element->ResultToPatch(values,nodesperelement,output_enum); + } + + /*Gather from all cpus*/ + ISSM_MPI_Allreduce((void*)values,(void*)allvalues,nbe*nodesperelement,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); + xDelete(values); + + if(save_results)results->AddResult(new GenericExternalResult(results->Size()+1,output_enum,allvalues,nbe,nodesperelement,step,time)); + xDelete(allvalues); } - Vector *vector_result = new Vector(size); - - /*Fill in vector*/ - for(int j=0;jSize();j++){ - Element* element=(Element*)elements->GetObjectByOffset(j); - element->ResultToVector(vector_result,output_enum); + else{ + + /*Allocate vector depending on interpolation*/ + switch(interpolation){ + case P0Enum: size = this->elements->NumberOfElements(); break; + case P1Enum: size = this->vertices->NumberOfVertices(); break; + default: _error_("Interpolation "< *vector_result = new Vector(size); + + /*Fill in vector*/ + for(int j=0;jSize();j++){ + Element* element=(Element*)elements->GetObjectByOffset(j); + element->ResultToVector(vector_result,output_enum); + } + vector_result->Assemble(); + + if(save_results)results->AddResult(new GenericExternalResult*>(results->Size()+1,output_enum,vector_result,step,time)); } - vector_result->Assemble(); - - if(save_results)results->AddResult(new GenericExternalResult*>(results->Size()+1,output_enum,vector_result,step,time)); - } - isvec = true; - break; - } + isvec = true; + break; + } + } + } Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 18871) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 18872) @@ -46,5 +46,6 @@ /*constructors, destructors: */ - FemModel(int argc,char** argv,ISSM_MPI_Comm comm_init); + FemModel(int argc,char** argv,ISSM_MPI_Comm comm_init,bool trace=false); + FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, ISSM_MPI_Comm incomm, int solution_type,IssmPDouble* X); ~FemModel(); @@ -52,5 +53,5 @@ void Echo(); FemModel* copy(); - void InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, const int solution_type); + void InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, const int solution_type,bool trace,IssmPDouble* X=NULL); void SolutionAnalysesList(int** panalyses,int* pnumanalyses,IoModel* iomodel,int solutiontype); void CleanUp(void);