Index: /issm/trunk/src/c/DataSet/DataSet.cpp =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.cpp (revision 961) +++ /issm/trunk/src/c/DataSet/DataSet.cpp (revision 962) @@ -173,4 +173,5 @@ if(object){ iterator = find(objects.begin(), objects.end(),object); + delete *iterator; objects.erase(iterator); } @@ -876,6 +877,6 @@ for ( object=objects.begin() ; object < objects.end(); object++ ){ delete (*object); - object=objects.erase(object); - } + } + objects.clear(); } Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 961) +++ /issm/trunk/src/c/Makefile.am (revision 962) @@ -32,4 +32,6 @@ ./objects/Friction.h\ ./objects/Friction.cpp\ + ./objects/DakotaPlugin.h\ + ./objects/DakotaPlugin.cpp\ ./objects/Node.h\ ./objects/Node.cpp\ @@ -266,4 +268,8 @@ ./ProcessParamsx/ProcessParamsx.cpp\ ./ProcessParamsx/ProcessParamsx.h\ + ./Qmux/Qmux.h\ + ./Qmux/Qmux.cpp\ + ./Qmux/SpawnCore.cpp\ + ./Qmux/SpawnCoreSerial.cpp\ ./FieldDepthAveragex/FieldDepthAveragex.cpp\ ./FieldDepthAveragex/FieldDepthAveragex.h\ @@ -536,4 +542,9 @@ ./ProcessParamsx/ProcessParamsx.cpp\ ./ProcessParamsx/ProcessParamsx.h\ + ./Qmux/Qmux.h\ + ./Qmux/Qmux.cpp\ + ./Qmux/DakotaResponses.cpp\ + ./Qmux/SpawnCore.cpp\ + ./Qmux/SpawnCoreParallel.cpp\ ./FieldDepthAveragex/FieldDepthAveragex.cpp\ ./FieldDepthAveragex/FieldDepthAveragex.h\ @@ -550,5 +561,4 @@ ./parallel/objectivefunctionC.cpp\ ./parallel/GradJCompute.cpp\ - ./parallel/SpawnCore.cpp\ ./parallel/ProcessResults.cpp\ ./parallel/prognostic_core.cpp\ @@ -556,6 +566,4 @@ ./parallel/transient_core_2d.cpp\ ./parallel/transient_core_3d.cpp\ - ./parallel/qmu.cpp\ - ./parallel/DakotaResponses.cpp\ ./parallel/OutputResults.cpp\ ./parallel/OutputControl.cpp Index: /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp (revision 961) +++ /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp (revision 962) @@ -21,4 +21,5 @@ /*create parameters common to all solutions: */ CreateParameters(pparameters,model,model_handle); + CreateParametersQmu(pparameters,model,model_handle); /*This is just a high level driver: */ @@ -39,5 +40,4 @@ CreateLoadsDiagnosticHoriz(ploads,model,model_handle); CreateParametersDiagnosticHoriz(pparameters,model,model_handle); - if(model->qmu_analysis)CreateParametersQmu(pparameters,model,model_handle); } @@ -78,5 +78,4 @@ CreateLoadsThermal(ploads,model,model_handle); CreateParametersThermal(pparameters,model,model_handle); - if(model->qmu_analysis)CreateParametersQmu(pparameters,model,model_handle); } @@ -94,5 +93,4 @@ CreateLoadsPrognostic(ploads,model,model_handle); CreateParametersPrognostic(pparameters,model,model_handle); - if(model->qmu_analysis)CreateParametersQmu(pparameters,model,model_handle); } @@ -100,4 +98,5 @@ throw ErrorException(__FUNCT__,exprintf("%s%s%s%s"," analysis_type: ",model->analysis_type," sub_analysis_type: ",model->sub_analysis_type," not supported yet!")); } + Index: /issm/trunk/src/c/ModelProcessorx/Model.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Model.cpp (revision 961) +++ /issm/trunk/src/c/ModelProcessorx/Model.cpp (revision 962) @@ -32,4 +32,5 @@ /*!initialize all pointers to 0: */ + model->name=NULL; model->repository=NULL; model->meshtype=NULL; @@ -40,5 +41,5 @@ model->numberofresponses=0; model->numberofvariables=0; - model->qmu_npart=0; + model->qmu_npart=0; model->numberofelements=0; model->numberofnodes=0; @@ -264,4 +265,5 @@ xfree((void**)&model->repository); xfree((void**)&model->meshtype); + xfree((void**)&model->name); xfree((void**)&model->analysis_type); xfree((void**)&model->sub_analysis_type); @@ -312,4 +314,5 @@ /*In ModelInit, we get all the data that is not difficult to get, and that is small: */ + ModelFetchData((void**)&model->name,NULL,NULL,model_handle,"name","String",NULL); ModelFetchData((void**)&model->analysis_type,NULL,NULL,model_handle,"analysis_type","String",NULL); ModelFetchData((void**)&model->sub_analysis_type,NULL,NULL,model_handle,"sub_analysis_type","String",NULL); Index: /issm/trunk/src/c/ModelProcessorx/Model.h =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Model.h (revision 961) +++ /issm/trunk/src/c/ModelProcessorx/Model.h (revision 962) @@ -13,4 +13,5 @@ struct Model { + char* name; char* repository; char* meshtype; Index: /issm/trunk/src/c/ModelProcessorx/Qmu/CreateParametersQmu.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Qmu/CreateParametersQmu.cpp (revision 961) +++ /issm/trunk/src/c/ModelProcessorx/Qmu/CreateParametersQmu.cpp (revision 962) @@ -11,4 +11,5 @@ #include "../../objects/objects.h" #include "../../shared/shared.h" +#include "../../include/macros.h" #include "../../MeshPartitionx/MeshPartitionx.h" #include "../Model.h" @@ -33,4 +34,14 @@ double* dakota_parameter=NULL; + //qmu files + char* qmuinname=NULL; + char* qmuerrname=NULL; + char* qmuoutname=NULL; + + #ifdef _SERIAL_ + mxArray* pfield=NULL; + mxArray* pfield2=NULL; + #endif + /*recover parameters : */ parameters=*pparameters; @@ -44,112 +55,173 @@ parameters->AddObject(param); - count++; - param= new Param(count,"qmu_npart",INTEGER); - param->SetInteger(model->qmu_npart); - parameters->AddObject(param); - - /*Deal with variables for qmu modeling: */ - variabledescriptors=(char**)xmalloc(model->numberofvariables*sizeof(char*)); - tag=(char*)xmalloc((strlen("variabledescriptori")+1)*sizeof(char)); - - /*Fetch descriptors: */ - for(i=0;inumberofvariables;i++){ - sprintf(tag,"%s%i","variabledescriptor",i); - ModelFetchData((void**)&descriptor,NULL,NULL,model_handle,tag,"String",NULL); - variabledescriptors[i]=descriptor; + + if(model->qmu_analysis){ + //name of qmu input, error and output files + qmuinname=(char*)xmalloc((strlen(model->name)+strlen(".qmu.in")+1)*sizeof(char)); + sprintf(qmuinname,"%s%s",model->name,".qmu.in"); + + count++; + param= new Param(count,"qmuinname",STRING); + param->SetString(qmuinname); + parameters->AddObject(param); + + qmuoutname=(char*)xmalloc((strlen(model->name)+strlen(".qmu.out")+1)*sizeof(char)); + sprintf(qmuoutname,"%s%s",model->name,".qmu.out"); + + count++; + param= new Param(count,"qmuoutname",STRING); + param->SetString(qmuoutname); + parameters->AddObject(param); + + qmuerrname=(char*)xmalloc((strlen(model->name)+strlen(".qmu.err")+1)*sizeof(char)); + sprintf(qmuerrname,"%s%s",model->name,".qmu.err"); + + count++; + param= new Param(count,"qmuerrname",STRING); + param->SetString(qmuerrname); + parameters->AddObject(param); + + //npart + count++; + param= new Param(count,"qmu_npart",INTEGER); + param->SetInteger(model->qmu_npart); + parameters->AddObject(param); + + /*Deal with variables for qmu modeling: */ + variabledescriptors=(char**)xmalloc(model->numberofvariables*sizeof(char*)); + + /*Fetch descriptors: logic varies if we are running parallel or serial. In parallel, qmumarshall + * took care of marshalling all the variable descriptors, so it's easy. In serial mode, + * the variables are in md.variables(md.ivar), as a strucuture: */ + + #ifdef _SERIAL_ + pfield=mxGetField(model_handle,0,"variabledescriptors"); + for(i=0;inumberofvariables;i++){ + pfield2=mxGetCell(pfield,i); + FetchData((void**)&descriptor,NULL,NULL,pfield2,"String",NULL); + variabledescriptors[i]=descriptor; + } + #else + tag=(char*)xmalloc((strlen("variabledescriptori")+1)*sizeof(char)); + for(i=0;inumberofvariables;i++){ + sprintf(tag,"%s%i","variabledescriptor",i); + ModelFetchData((void**)&descriptor,NULL,NULL,model_handle,tag,"String",NULL); + variabledescriptors[i]=descriptor; + } + #endif + + /*Ok, we have all the variable descriptors. Build a parameter with it: */ + count++; + param= new Param(count,"variabledescriptors",STRINGARRAY); + param->SetStringArray(variabledescriptors,model->numberofvariables); + parameters->AddObject(param); + + + /*Deal with responses and partition for qmu modeling: */ + responsedescriptors=(char**)xmalloc(model->numberofresponses*sizeof(char*)); + + /*Fetch descriptors: */ + #ifdef _SERIAL_ + pfield=mxGetField(model_handle,0,"responsedescriptors"); + for(i=0;inumberofresponses;i++){ + pfield2=mxGetCell(pfield,i); + FetchData((void**)&descriptor,NULL,NULL,pfield2,"String",NULL); + responsedescriptors[i]=descriptor; + } + #else + xfree((void**)&tag); + tag=(char*)xmalloc((strlen("responsedescriptori")+1)*sizeof(char)); + + for(i=0;inumberofresponses;i++){ + sprintf(tag,"%s%i","responsedescriptor",i); + ModelFetchData((void**)&descriptor,NULL,NULL,model_handle,tag,"String",NULL); + responsedescriptors[i]=descriptor; + } + #endif + + /*Ok, we have all the response descriptors. Build a parameter with it: */ + count++; + param= new Param(count,"responsedescriptors",STRINGARRAY); + param->SetStringArray(responsedescriptors,model->numberofresponses); + parameters->AddObject(param); + + #ifdef _DEBUG_ + for(i=0;inumberofvariables;i++){ + _printf_("variable descriptor %s\n",variabledescriptors[i]); + } + + for(i=0;inumberofresponses;i++){ + _printf_("response descriptor %s\n",responsedescriptors[i]); + } + #endif + + /*partition grids in model->qmu_npart parts: */ + + if(strcmp(model->meshtype,"2d")==0){ + ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat"); + elements_width=3; //tria elements + } + else{ + ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat"); + elements_width=6; //penta elements + } + + MeshPartitionx(&epart, &part,model->numberofelements,model->numberofnodes,model->elements, model->numberofelements2d,model->numberofnodes2d,model->elements2d,model->numlayers,elements_width, model->meshtype,model->qmu_npart); + + dpart=(double*)xmalloc(model->numberofnodes*sizeof(double)); + for(i=0;inumberofnodes;i++)dpart[i]=part[i]; + + count++; + param= new Param(count,"qmu_part",DOUBLEVEC); + param->SetDoubleVec(dpart,model->numberofnodes,1); + parameters->AddObject(param); + + /*Ok, now if any of the variables input from Dakota are distributed, we are going to need the parameters: */ + for(i=0;inumberofvariables;i++){ + + descriptor=variabledescriptors[i]; + + if ((strcmp(descriptor,"thickness")==0) || + (strcmp(descriptor,"drag") ==0) + ){ + + //Fetch data: + ModelFetchData((void**)&dakota_parameter,NULL,NULL,model_handle,descriptor,"Matrix","Mat"); + + //Add parameter + count++; + param= new Param(count,descriptor,DOUBLEVEC); + param->SetDoubleVec(dakota_parameter,model->numberofnodes,1); + parameters->AddObject(param); + + //free data + xfree((void**)&dakota_parameter); + + } + } + + /*Free data: */ + xfree((void**)&tag); + for(i=0;inumberofresponses;i++){ + char* descriptor=responsedescriptors[i]; + xfree((void**)&descriptor); + } + xfree((void**)&responsedescriptors); + + for(i=0;inumberofvariables;i++){ + char* descriptor=variabledescriptors[i]; + xfree((void**)&descriptor); + } + xfree((void**)&variabledescriptors); + + xfree((void**)&model->elements); + xfree((void**)&model->elements2d); + xfree((void**)&epart); + xfree((void**)&part); + xfree((void**)&dpart); + xfree((void**)&qmuinname); + xfree((void**)&qmuerrname); + xfree((void**)&qmuoutname); } - - /*Ok, we have all the variable descriptors. Build a parameter with it: */ - count++; - param= new Param(count,"variabledescriptors",STRINGARRAY); - param->SetStringArray(variabledescriptors,model->numberofvariables); - parameters->AddObject(param); - - xfree((void**)&tag); - - /*Deal with responses and partition for qmu modeling: */ - responsedescriptors=(char**)xmalloc(model->numberofresponses*sizeof(char*)); - tag=(char*)xmalloc((strlen("responsedescriptori")+1)*sizeof(char)); - - /*Fetch descriptors: */ - for(i=0;inumberofresponses;i++){ - sprintf(tag,"%s%i","responsedescriptor",i); - ModelFetchData((void**)&descriptor,NULL,NULL,model_handle,tag,"String",NULL); - responsedescriptors[i]=descriptor; - } - - /*Ok, we have all the response descriptors. Build a parameter with it: */ - count++; - param= new Param(count,"responsedescriptors",STRINGARRAY); - param->SetStringArray(responsedescriptors,model->numberofresponses); - parameters->AddObject(param); - - - #ifdef _PARALLEL_ - /*partition grids in model->qmu_npart parts: */ - - if(strcmp(model->meshtype,"2d")==0){ - ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat"); - elements_width=3; //tria elements - } - else{ - ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat"); - elements_width=6; //penta elements - } - - MeshPartitionx(&epart, &part,model->numberofelements,model->numberofnodes,model->elements, model->numberofelements2d,model->numberofnodes2d,model->elements2d,model->numlayers,elements_width, model->meshtype,model->qmu_npart); - - dpart=(double*)xmalloc(model->numberofnodes*sizeof(double)); - for(i=0;inumberofnodes;i++)dpart[i]=part[i]; - - count++; - param= new Param(count,"qmu_part",DOUBLEVEC); - param->SetDoubleVec(dpart,model->numberofnodes,1); - parameters->AddObject(param); - #endif - - /*Ok, now if any of the variables input from Dakota are distributed, we are going to need the parameters: */ - for(i=0;inumberofvariables;i++){ - - descriptor=variabledescriptors[i]; - - if ((strcmp(descriptor,"thickness")==0) || - (strcmp(descriptor,"drag") ==0) - ){ - - //Fetch data: - ModelFetchData((void**)&dakota_parameter,NULL,NULL,model_handle,descriptor,"Matrix","Mat"); - - //Add parameter - count++; - param= new Param(count,descriptor,DOUBLEVEC); - param->SetDoubleVec(dakota_parameter,model->numberofnodes,1); - parameters->AddObject(param); - - //free data - xfree((void**)&dakota_parameter); - - } - } - - /*Free data: */ - xfree((void**)&tag); - for(i=0;inumberofresponses;i++){ - char* descriptor=responsedescriptors[i]; - xfree((void**)&descriptor); - } - xfree((void**)&responsedescriptors); - - for(i=0;inumberofvariables;i++){ - char* descriptor=variabledescriptors[i]; - xfree((void**)&descriptor); - } - xfree((void**)&variabledescriptors); - - xfree((void**)&model->elements); - xfree((void**)&model->elements2d); - xfree((void**)&epart); - xfree((void**)&part); - xfree((void**)&dpart); Index: /issm/trunk/src/c/Qmux/DakotaResponses.cpp =================================================================== --- /issm/trunk/src/c/Qmux/DakotaResponses.cpp (revision 962) +++ /issm/trunk/src/c/Qmux/DakotaResponses.cpp (revision 962) @@ -0,0 +1,190 @@ +/*!\file: DakotaResponses.cpp + * \brief compute dakota responses, using a list of response descriptors. + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../DataSet/DataSet.h" +#include "../shared/shared.h" + +#undef __FUNCT__ +#define __FUNCT__ "DakotaResponses" +void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,FemModel* femmodels,DataSet* results,int analysis_type,int sub_analysis_type){ + + int i,j; + int found=0; + char* response_descriptor=NULL; + int numberofnodes; + FemModel* fem=NULL; + extern int my_rank; + + /*recover first model: */ + fem=femmodels; + + /*some data needed across the responses: */ + found=fem->parameters->FindParam((void*)&numberofnodes,"numberofnodes"); + if(!found)throw ErrorException(__FUNCT__," could not find numberofnodes in fem model"); + + + for(i=0;iFindResult((void*)&vel,"vel"); + if(!found)throw ErrorException(__FUNCT__," could not find vel to compute min_vel"); + + min_vel=vel[0]; + for(j=1;jFindResult((void*)&vel,"vel"); + if(!found)throw ErrorException(__FUNCT__," could not find vel to compute max_vel"); + + max_vel=vel[0]; + for(j=1;jmax_vel)max_vel=vel[j]; + } + if(my_rank==0)responses[i]=max_vel; + } + else if(strcmp(response_descriptor,"min_vx")==0){ + double* vx=NULL; + double min_vx=0; + + found=results->FindResult((void*)&vx,"vx"); + if(!found)throw ErrorException(__FUNCT__," could not find vx to compute min_vx"); + + min_vx=vx[0]; + for(j=1;jFindResult((void*)&vx,"vx"); + if(!found)throw ErrorException(__FUNCT__," could not find vx to compute max_vx"); + + max_vx=vx[0]; + for(j=1;jmax_vx)max_vx=vx[j]; + } + if(my_rank==0)responses[i]=max_vx; + } + else if(strcmp(response_descriptor,"max_abs_vx")==0){ + double* vx=NULL; + double max_abs_vx=0; + + found=results->FindResult((void*)&vx,"vx"); + if(!found)throw ErrorException(__FUNCT__," could not find vx to compute max_abs_vx"); + + max_abs_vx=fabs(vx[0]); + for(j=1;jmax_abs_vx)max_abs_vx=fabs(vx[j]); + } + if(my_rank==0)responses[i]=max_abs_vx; + } + else if(strcmp(response_descriptor,"min_vy")==0){ + double* vy=NULL; + double min_vy=0; + + found=results->FindResult((void*)&vy,"vy"); + if(!found)throw ErrorException(__FUNCT__," could not find vy to compute min_vy"); + + min_vy=vy[0]; + for(j=1;jFindResult((void*)&vy,"vy"); + if(!found)throw ErrorException(__FUNCT__," could not find vy to compute max_vy"); + + max_vy=vy[0]; + for(j=1;jmax_vy)max_vy=vy[j]; + } + if(my_rank==0)responses[i]=max_vy; + } + else if(strcmp(response_descriptor,"max_abs_vy")==0){ + double* vy=NULL; + double max_abs_vy=0; + + found=results->FindResult((void*)&vy,"vy"); + if(!found)throw ErrorException(__FUNCT__," could not find vy to compute max_abs_vy"); + + max_abs_vy=fabs(vy[0]); + for(j=1;jmax_abs_vy)max_abs_vy=fabs(vy[j]); + } + if(my_rank==0)responses[i]=max_abs_vy; + } + else if(strcmp(response_descriptor,"min_vz")==0){ + double* vz=NULL; + double min_vz=0; + + found=results->FindResult((void*)&vz,"vz"); + if(!found)throw ErrorException(__FUNCT__," could not find vz to compute min_vz"); + + min_vz=vz[0]; + for(j=1;jFindResult((void*)&vz,"vz"); + if(!found)throw ErrorException(__FUNCT__," could not find vz to compute max_vz"); + + max_vz=vz[0]; + for(j=1;jmax_vz)max_vz=vz[j]; + } + if(my_rank==0)responses[i]=max_vz; + } + else if(strcmp(response_descriptor,"max_abs_vz")==0){ + double* vz=NULL; + double max_abs_vz=0; + + found=results->FindResult((void*)&vz,"vz"); + if(!found)throw ErrorException(__FUNCT__," could not find vz to compute max_abs_vz"); + + max_abs_vz=fabs(vz[0]); + for(j=1;jmax_abs_vz)max_abs_vz=fabs(vz[j]); + } + if(my_rank==0)responses[i]=max_abs_vz; + } + else{ + if(my_rank==0)printf("%s%s%s"," response descriptor : ",response_descriptor," not supported yet!"); + throw ErrorException(__FUNCT__,exprintf("%s%s%s"," response descriptor : ",response_descriptor," not supported yet!")); + } + } + +} Index: /issm/trunk/src/c/Qmux/Qmux.cpp =================================================================== --- /issm/trunk/src/c/Qmux/Qmux.cpp (revision 962) +++ /issm/trunk/src/c/Qmux/Qmux.cpp (revision 962) @@ -0,0 +1,162 @@ +/*!\file: Qmux.cpp + * \brief: wrapper to the Dakota capabilities. qmu fires up Dakota, and registers a Dakota Pluggin + * which will be in charge of running the solution sequences repeteadly, to garner statistics. + * + * This routine deals with running ISSM and Dakota in library mode. In library mode, Dakota does not + * run as an execuatble. Its capabilities are linked into the ISSM software. ISSM calls dakota routines + * directly from the dakota library. qmu.cpp is the code that is in charge of calling those routines. + * + * Dakota has its own way of running in parallel (for embarassingly parallel jobs). We do not want that, + * as ISSM knows exactly how to run "really parallel" jobs that use all CPUS. To bypass Dakota's parallelism, + * we overloaded the constructor for the parallel library (see the Dakota patch in the externalpackages/dakota + * directory). This overloaded constructor fires up Dakota serially on CPU 0 only! We take care of broadcasting + * to the other CPUS, hence ISSM is running in parallel, and Dakota serially on CPU0. + * + * Now, how does CPU 0 drive all other CPUS to carry out sensitivity analysese? By synchronizing its call to + * our ISSM cores (diagnostic_core, thermal_core, transient_core, etc ...) on CPU 0 with all other CPUS. + * This explains the structure of qmu.cpp, where cpu 0 runs Dakota, the Dakota pluggin fires up SpawnCore.cpp, + * while the other CPUS are waiting for a broadcast from CPU0, once they get it, they also fire up + * SpawnCore. In the end, SpawnCore is fired up on all CPUS, with CPU0 having Dakota inputs, that it will + * broacast to other CPUS. + * + * Now, how does dakota call the SpawnCore routine? The SpawnCore is embedded into the DakotaPlugin object + * which is derived from the Direct Interface Dakota objct. This is the only way to run Dakota in library + * mode (see their developper guide for more info). Dakota registers the DakotaPlugin object into its own + * database, and calls on the embedded SpawnCore from CPU0. + * + */ + + + +#ifdef HAVE_CONFIG_H +#include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + + +#undef __FUNCT__ +#define __FUNCT__ "Qmux" + +#include "./Qmux.h" + +#include "../shared/shared.h" +#include "../include/macros.h" +#include "../toolkits/toolkits.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../objects/DakotaPlugin.h" + +#ifdef _HAVE_DAKOTA_ //only works if dakota library has been compiled in. +#include "ParallelLibrary.H" +#include "ProblemDescDB.H" +#include "DakotaStrategy.H" +#include "DakotaModel.H" +#include "DakotaInterface.H" + +#endif + +#ifdef _SERIAL_ +void Qmux(mxArray* femmodels,mxArray* inputs,int analysis_type,int sub_analysis_type,char* dakota_input_file,char* dakota_output_file,char* dakota_error_file){ +#else +void Qmux(FemModel* femmodels,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ +#endif + + + #ifdef _HAVE_DAKOTA_ //only works if dakota library has been compiled in. + + #ifdef _PARALLEL_ + char* dakota_input_file=NULL; + char* dakota_output_file=NULL; + char* dakota_error_file=NULL; + extern int my_rank; + #endif + int status=0; + Dakota::ModelLIter ml_iter; + + #ifdef _PARALLEL_ + /*Recover dakota_input_file, dakota_output_file and dakota_error_file, in the parameters dataset in parallel */ + femmodels->parameters->FindParam((void*)&dakota_input_file,"qmuinname"); + femmodels->parameters->FindParam((void*)&dakota_output_file,"qmuoutname"); + femmodels->parameters->FindParam((void*)&dakota_error_file,"qmuerrname"); + #endif + + #ifdef _PARALLEL_ + if(my_rank==0){ + #endif + + // Instantiate/initialize the parallel library and problem description + // database objects. + Dakota::ParallelLibrary parallel_lib("serial"); + Dakota::ProblemDescDB problem_db(parallel_lib); + + // Manage input file parsing, output redirection, and restart processing + // without a CommandLineHandler. This version relies on parsing of an + // input file. + problem_db.manage_inputs(dakota_input_file); + // specify_outputs_restart() is only necessary if specifying non-defaults + parallel_lib.specify_outputs_restart(dakota_output_file,dakota_error_file,NULL,NULL); + + // Instantiate the Strategy object (which instantiates all Model and + // Iterator objects) using the parsed information in problem_db. + Dakota::Strategy selected_strategy(problem_db); + + // convenience function for iterating over models and performing any + // interface plug-ins + Dakota::ModelList& models = problem_db.model_list(); + + for (ml_iter = models.begin(); ml_iter != models.end(); ml_iter++) { + + Dakota::Interface& interface = ml_iter->interface(); + + //set DB nodes to the existing Model specification + problem_db.set_db_model_nodes(ml_iter->model_id()); + + // Serial case: plug in derived Interface object without an analysisComm + interface.assign_rep(new SIM::DakotaPlugin(problem_db,(void*)femmodels,(void*)inputs,analysis_type,sub_analysis_type), false); + } + + // Execute the strategy + problem_db.lock(); // prevent run-time DB queries + selected_strategy.run_strategy(); + + #ifdef _PARALLEL_ + //Warn other cpus that we are done running the dakota iterator + MPI_Bcast(&status,1,MPI_INT,0,MPI_COMM_WORLD); + #endif + #ifdef _PARALLEL_ + } + else{ + + for(;;){ + + SpawnCore(NULL,0, NULL,NULL,0,femmodels,inputs,analysis_type,sub_analysis_type,0); + + //Figure out if cpu 0 is done iterating + MPI_Bcast(&status,1,MPI_INT,0,MPI_COMM_WORLD); + + //yes! + if(status){ + break; //yes, we are done + } + else{ + /*We are not done: */ + + + } + } + + + } + #endif + + /*Free ressources:*/ + #ifdef _PARALLEL_ + xfree((void**)&dakota_input_file); + xfree((void**)&dakota_error_file); + xfree((void**)&dakota_output_file); + #endif + + + #endif + +} Index: /issm/trunk/src/c/Qmux/Qmux.h =================================================================== --- /issm/trunk/src/c/Qmux/Qmux.h (revision 962) +++ /issm/trunk/src/c/Qmux/Qmux.h (revision 962) @@ -0,0 +1,23 @@ +/*!\file: Qmux.h + * \brief header file for Qmu engine using Dakota + */ + +#ifndef _QMUX_H +#define _QMUX_H + +#include "../DataSet/DataSet.h" + +/* local prototypes: */ +void SpawnCore(double* responses, int numresponses, double* variables, char** variables_descriptors,int numvariables, void* femmodels,void* inputs,int analysis_type,int sub_analysis_type,int counter); +#ifdef _SERIAL_ +void Qmux(mxArray* femmodels,mxArray* inputs,int analysis_type,int sub_analysis_type,char* dakota_input_file,char* dakota_output_file,char* dakota_error_file); +void SpawnCoreSerial(double* responses, int numresponses, double* variables, char** variables_descriptors,int numvariables, mxArray* femmodels,mxArray* inputs,int analysis_type,int sub_analysis_type,int counter); +#else +void Qmux(FemModel* femmodels,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); +void SpawnCoreParallel(double* responses, int numresponses, double* variables, char** variables_descriptors,int numvariables, FemModel* femmodels,ParameterInputs* inputs,int analysis_type,int sub_analysis_type,int counter); +void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,FemModel* femmodels, DataSet* results,int analysis_type,int sub_analysis_type); +#endif + + +#endif /* _QMUX_H */ + Index: /issm/trunk/src/c/Qmux/SpawnCore.cpp =================================================================== --- /issm/trunk/src/c/Qmux/SpawnCore.cpp (revision 962) +++ /issm/trunk/src/c/Qmux/SpawnCore.cpp (revision 962) @@ -0,0 +1,32 @@ +/*!\file: SpawnCore.cpp + * \brief: branch into SpawnCoreMatlab and SpawnCoreParallel. + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#undef __FUNCT__ +#define __FUNCT__ "SpawnCore" + +#include "../objects/objects.h" +#include "../io/io.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../shared/shared.h" +#include "./Qmux.h" +#include "../include/macros.h" + +void SpawnCore(double* responses, int numresponses, double* variables, char** variables_descriptors,int numvariables, void* femmodels,void* inputs,int analysis_type,int sub_analysis_type,int counter){ + + + /*Branch into a serial SpawnCore and a parallel SpawnCore: */ + #ifdef _SERIAL_ + SpawnCoreSerial(responses, numresponses, variables, variables_descriptors,numvariables, (mxArray*)femmodels, (mxArray*)inputs,analysis_type,sub_analysis_type,counter); + #else + SpawnCoreParallel(responses, numresponses, variables, variables_descriptors,numvariables, (FemModel*)femmodels, (ParameterInputs*)inputs,analysis_type,sub_analysis_type,counter); + #endif + + +} Index: /issm/trunk/src/c/Qmux/SpawnCoreParallel.cpp =================================================================== --- /issm/trunk/src/c/Qmux/SpawnCoreParallel.cpp (revision 962) +++ /issm/trunk/src/c/Qmux/SpawnCoreParallel.cpp (revision 962) @@ -0,0 +1,180 @@ +/*!\file: SpawnCoreParallel.cpp + * \brief: run core ISSM solution using Dakota inputs coming from CPU 0. + * \sa qmu.cpp DakotaPlugin.cpp + * + * This routine needs to be understood simultaneously with qmu.cpp and DakotaPlugin. + * SpawnCoreParallel is called by all CPUS, with CPU 0 holding Dakota variable values, along + * with variable descriptors. + * + * SpawnCoreParallel takes care of broadcasting the variables and their descriptors across the MPI + * ring. Once this is done, we use the variables to modify the inputs for the solution core. + * For ex, if "rho_ice" is provided, for ex 920, we include "rho_ice" in the inputs, then + * call the core with the modified inputs. This is the way we get Dakota to explore the parameter + * spce of the core. + * + * Once the core is called, we process the results of the core, and using the processed results, + * we compute response functions. The responses are computed on all CPUS, but they are targeted + * for CPU 0, which will get these values back to the Dakota engine. + * + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#undef __FUNCT__ +#define __FUNCT__ "SpawnCoreParallel" + +#include "../objects/objects.h" +#include "../io/io.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../shared/shared.h" +#include "./Qmux.h" +#include "../include/macros.h" +#include "../parallel/parallel.h" + +void SpawnCoreParallel(double* responses, int numresponses, double* variables, char** variables_descriptors,int numvariables, FemModel* femmodels,ParameterInputs* inputs,int analysis_type,int sub_analysis_type,int counter){ + + int i; + + /*output from core solutions: */ + DataSet* results=NULL; + + char** responses_descriptors=NULL; + Param* param=NULL; + char* string=NULL; + int string_length; + double* qmu_part=NULL; + int qmu_npart; + int debug=0; + + extern int my_rank; + + + /*synchronize all cpus, as CPU 0 is probably late (it is starting the entire dakota strategy!) : */ + MPI_Barrier(MPI_COMM_WORLD); + + /*some parameters needed: */ + femmodels[0].parameters->FindParam((void*)&debug,"debug"); + + /*First off, recover the response descriptors for the response functions: */ + param=(Param*)femmodels[0].parameters->FindParamObject("responsedescriptors"); + if(!param)throw ErrorException(__FUNCT__," could not find response descriptors!"); + + param->GetParameterValue(&responses_descriptors); + + /*Recover partitioning for dakota: */ + femmodels[0].parameters->FindParam((void*)&qmu_npart,"qmu_npart"); + femmodels[0].parameters->FindParam((void*)&qmu_part,"qmu_part"); + #ifdef _DEBUG_ + for(i=0;iUpdateFromDakota(variables,variables_descriptors,numvariables,femmodels[0].parameters,qmu_part,qmu_npart); //femmodel #0 is the one holding the parameters for Dakota. + + /*Run the analysis core solution sequence, with the updated inputs: */ + if(analysis_type==DiagnosticAnalysisEnum()){ + + if(debug)_printf_("Starting diagnostic core\n"); + + diagnostic_core(results,femmodels,inputs); + + } + else if(analysis_type==ThermalAnalysisEnum()){ + + if(debug)_printf_("Starting thermal core\n"); + thermal_core(results,femmodels,inputs); + + } + else if(analysis_type==PrognosticAnalysisEnum()){ + + if(debug)_printf_("Starting prognostic core\n"); + prognostic_core(results,femmodels,inputs); + + } + else if(analysis_type==TransientAnalysisEnum()){ + + if(debug)_printf_("Starting transient core\n"); + transient_core(results,femmodels,inputs); + + } + else throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s"," analysis_type ",analysis_type," and sub_analysis_type ",sub_analysis_type," not supported yet!")); + + + + /*Now process the outputs, before computing the dakota responses: */ + if(debug)_printf_("process results:\n"); + ProcessResults(&results,femmodels,analysis_type); + + + /*compute responses on cpu 0: dummy for now! */ + if(debug)_printf_("compute dakota responses:\n"); + DakotaResponses(responses,responses_descriptors,numresponses,femmodels,results,analysis_type,sub_analysis_type); + + /*Free ressources:*/ + delete results; + + //variables only on cpu != 0 + if(my_rank!=0){ + xfree((void**)&variables); + for(i=0;iGetParameterValue((void*)&string); mxSetField( dataref, 0, param->GetParameterName(),mxCreateString(string)); + break; + + case STRINGARRAY: + param->GetParameterValue((void*)&stringarray); + M=param->GetM(); + dims[0]=M; + dims[1]=1; + pfield=mxCreateCellArray(2,dims); + for(k=0;kGetParameterName(),pfield); break; Index: /issm/trunk/src/c/issm.h =================================================================== --- /issm/trunk/src/c/issm.h (revision 961) +++ /issm/trunk/src/c/issm.h (revision 962) @@ -50,4 +50,5 @@ #include "./ComputePressurex/ComputePressurex.h" #include "./FieldExtrudex/FieldExtrudex.h" +#include "./Qmux/Qmux.h" Index: /issm/trunk/src/c/objects/DakotaPlugin.cpp =================================================================== --- /issm/trunk/src/c/objects/DakotaPlugin.cpp (revision 961) +++ /issm/trunk/src/c/objects/DakotaPlugin.cpp (revision 962) @@ -26,5 +26,5 @@ #include "../include/macros.h" #include "../objects/objects.h" -#include "../parallel/parallel.h" +#include "../Qmux/Qmux.h" /*Standard includes: */ @@ -44,5 +44,6 @@ //constructor -DakotaPlugin::DakotaPlugin(const Dakota::ProblemDescDB& problem_db,FemModel* in_femmodels, ParameterInputs* in_inputs, int in_analysis_type, int in_sub_analysis_type):Dakota::DirectApplicInterface(problem_db){ +DakotaPlugin::DakotaPlugin(const Dakota::ProblemDescDB& problem_db,void* in_femmodels, void* in_inputs, int in_analysis_type, int in_sub_analysis_type):Dakota::DirectApplicInterface(problem_db){ + femmodels=in_femmodels; @@ -90,5 +91,5 @@ /*run core solution: */ - SpawnCore(responses,variables,variable_descriptors,numACV,femmodels,inputs,analysis_type,sub_analysis_type,counter); + SpawnCore(responses,numFns, variables,variable_descriptors,numACV,femmodels,inputs,analysis_type,sub_analysis_type,counter); /*populate responses: */ @@ -97,6 +98,8 @@ } + #ifdef _PARALLEL_ /*warn other cpus that we are done running this iteration: */ MPI_Bcast(&status,1,MPI_INT,0,MPI_COMM_WORLD); + #endif /*Free ressources:*/ Index: /issm/trunk/src/c/objects/DakotaPlugin.h =================================================================== --- /issm/trunk/src/c/objects/DakotaPlugin.h (revision 961) +++ /issm/trunk/src/c/objects/DakotaPlugin.h (revision 962) @@ -18,13 +18,14 @@ public: - DakotaPlugin(const Dakota::ProblemDescDB& problem_db,FemModel* femmodels, ParameterInputs* inputs, int analysis_type, int sub_analysis_type); - ~DakotaPlugin(); - - /*these fields are use by core solutions: */ - FemModel* femmodels; - ParameterInputs* inputs; - int analysis_type; - int sub_analysis_type; - int counter; + DakotaPlugin(const Dakota::ProblemDescDB& problem_db,void* femmodels, void* inputs, int analysis_type, int sub_analysis_type); + ~DakotaPlugin(); + + /*these fields are used by core solutions: */ + void* femmodels; + void* inputs; + + int analysis_type; + int sub_analysis_type; + int counter; protected: Index: /issm/trunk/src/c/objects/Param.cpp =================================================================== --- /issm/trunk/src/c/objects/Param.cpp (revision 961) +++ /issm/trunk/src/c/objects/Param.cpp (revision 962) @@ -58,5 +58,4 @@ case STRINGARRAY: - for(i=0;iAddObject(input); @@ -329,4 +329,14 @@ } + #ifdef _DEBUG_ + PetscSynchronizedPrintf(MPI_COMM_WORLD,"Parameter vetor:"); + PetscSynchronizedFlush(MPI_COMM_WORLD); + for(k=0;kAdd(root,parameter,1,numberofnodes); Index: /issm/trunk/src/c/parallel/OutputControl.cpp =================================================================== --- /issm/trunk/src/c/parallel/OutputControl.cpp (revision 961) +++ /issm/trunk/src/c/parallel/OutputControl.cpp (revision 962) @@ -6,4 +6,5 @@ #include "../EnumDefinitions/EnumDefinitions.h" #include "../objects/objects.h" +#include "../shared/shared.h" #ifdef HAVE_CONFIG_H @@ -24,5 +25,5 @@ int gsize; - /*recover numberofnodes: */ + /*recover parameters: */ fem->parameters->FindParam((void*)&numberofnodes,"numberofnodes"); gsize=fem->nodesets->GetGSize(); Index: sm/trunk/src/c/parallel/SpawnCore.cpp =================================================================== --- /issm/trunk/src/c/parallel/SpawnCore.cpp (revision 961) +++ (revision ) @@ -1,183 +1,0 @@ -/*!\file: SpawnCore.cpp - * \brief: run core ISSM solution using Dakota inputs coming from CPU 0. - * \sa qmu.cpp DakotaPlugin.cpp - * - * This routine needs to be understood simultaneously with qmu.cpp and DakotaPlugin. - * SpawnCore is called by all CPUS, with CPU 0 holding Dakota variable values, along - * with variable descriptors. - * - * SpawnCore takes care of broadcasting the variables and their descriptors across the MPI - * ring. Once this is done, we use the variables to modify the inputs for the solution core. - * For ex, if "rho_ice" is provided, for ex 920, we include "rho_ice" in the inputs, then - * call the core with the modified inputs. This is the way we get Dakota to explore the parameter - * spce of the core. - * - * Once the core is called, we process the results of the core, and using the processed results, - * we compute response functions. The responses are computed on all CPUS, but they are targeted - * for CPU 0, which will get these values back to the Dakota engine. - * - */ - -#ifdef HAVE_CONFIG_H - #include "config.h" -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif - -#undef __FUNCT__ -#define __FUNCT__ "SpawnCore" - -#include "../objects/objects.h" -#include "../io/io.h" -#include "../EnumDefinitions/EnumDefinitions.h" -#include "../shared/shared.h" -#include "./parallel.h" -#include "../include/macros.h" - -void SpawnCore(double* responses, double* variables, char** variables_descriptors,int numvariables, FemModel* femmodels,ParameterInputs* inputs,int analysis_type,int sub_analysis_type,int counter){ - - int i; - - /*output from core solutions: */ - DataSet* results=NULL; - - char** responses_descriptors=NULL; - int numresponses; - Param* param=NULL; - char* string=NULL; - int string_length; - double* qmu_part=NULL; - int qmu_npart; - int debug=0; - - extern int my_rank; - - - /*synchronize all cpus, as CPU 0 is probably late (it is starting the entire dakota strategy!) : */ - MPI_Barrier(MPI_COMM_WORLD); - - /*some parameters needed: */ - femmodels[0].parameters->FindParam((void*)&debug,"debug"); - - /*First off, recover the response descriptors for the response functions: */ - param=(Param*)femmodels[0].parameters->FindParamObject("responsedescriptors"); - if(!param)throw ErrorException(__FUNCT__," could not find response descriptors!"); - - numresponses=param->GetM(); - param->GetParameterValue(&responses_descriptors); - - /*Recover partitioning for dakota: */ - femmodels[0].parameters->FindParam((void*)&qmu_npart,"qmu_npart"); - femmodels[0].parameters->FindParam((void*)&qmu_part,"qmu_part"); - #ifdef _DEBUG_ - for(i=0;iUpdateFromDakota(variables,variables_descriptors,numvariables,femmodels[0].parameters,qmu_part,qmu_npart); //femmodel #0 is the one holding the parameters for Dakota. - - /*Run the analysis core solution sequence, with the updated inputs: */ - if(analysis_type==DiagnosticAnalysisEnum()){ - - if(debug)_printf_("Starting diagnostic core\n"); - - diagnostic_core(results,femmodels,inputs); - - } - else if(analysis_type==ThermalAnalysisEnum()){ - - if(debug)_printf_("Starting thermal core\n"); - thermal_core(results,femmodels,inputs); - - } - else if(analysis_type==PrognosticAnalysisEnum()){ - - if(debug)_printf_("Starting prognostic core\n"); - prognostic_core(results,femmodels,inputs); - - } - else if(analysis_type==TransientAnalysisEnum()){ - - if(debug)_printf_("Starting transient core\n"); - transient_core(results,femmodels,inputs); - - } - else throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s"," analysis_type ",analysis_type," and sub_analysis_type ",sub_analysis_type," not supported yet!")); - - - - /*Now process the outputs, before computing the dakota responses: */ - if(debug)_printf_("process results:\n"); - ProcessResults(&results,femmodels,analysis_type); - - - /*compute responses on cpu 0: dummy for now! */ - if(debug)_printf_("compute dakota responses:\n"); - DakotaResponses(responses,responses_descriptors,numresponses,femmodels,results,analysis_type,sub_analysis_type); - - /*Free ressources:*/ - delete results; - - //variables only on cpu != 0 - if(my_rank!=0){ - xfree((void**)&variables); - for(i=0;iAdd(control_type,param_g,2,numberofnodes); @@ -161,5 +161,5 @@ OutputControl(u_g,param_g,J,nsteps,&femmodel,outputfilename); _printf_("%s\n"," done."); - } + }*/ _printf_("%s%i%s%g\n"," value of misfit J after optimization #",n,": ",J[n]); Index: /issm/trunk/src/c/parallel/diagnostic.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic.cpp (revision 961) +++ /issm/trunk/src/c/parallel/diagnostic.cpp (revision 962) @@ -23,7 +23,4 @@ char* outputfilename=NULL; char* lockname=NULL; - char* qmuinname=NULL; - char* qmuoutname=NULL; - char* qmuerrname=NULL; int numberofnodes; int qmu_analysis=0; @@ -57,7 +54,4 @@ outputfilename=argv[3]; lockname=argv[4]; - qmuinname=argv[5]; - qmuoutname=argv[6]; - qmuerrname=argv[7]; /*Open handle to data on disk: */ @@ -105,5 +99,5 @@ #ifdef _HAVE_DAKOTA_ - qmu(qmuinname,qmuoutname,qmuerrname,&femmodels[0],inputs,DiagnosticAnalysisEnum(),NoneAnalysisEnum()); + Qmux(&femmodels[0],inputs,DiagnosticAnalysisEnum(),NoneAnalysisEnum()); #else throw ErrorException(__FUNCT__," Dakota not present, cannot do qmu!"); Index: /issm/trunk/src/c/parallel/parallel.h =================================================================== --- /issm/trunk/src/c/parallel/parallel.h (revision 961) +++ /issm/trunk/src/c/parallel/parallel.h (revision 962) @@ -42,7 +42,4 @@ void CreateFemModel(FemModel* femmodel,ConstDataHandle MODEL,char* analysis_type,char* sub_analysis_type); //int BatchDebug(Mat* Kgg,Vec* pg,FemModel* femmodel,char* filename); -void qmu(const char* dakota_input_file,const char* dakota_output_file,const char* dakota_error_file,FemModel* femmodels,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); -void SpawnCore(double* responses,double* variables,char** variable_descriptors,int numvariables, FemModel* femmodels,ParameterInputs* inputs,int analysis_type,int sub_analysis_type,int counter); -void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,FemModel* femmodels, DataSet* results,int analysis_type,int sub_analysis_type); void ProcessResults(DataSet** presults,FemModel* fems,int analysis_type); Index: /issm/trunk/src/c/parallel/prognostic.cpp =================================================================== --- /issm/trunk/src/c/parallel/prognostic.cpp (revision 961) +++ /issm/trunk/src/c/parallel/prognostic.cpp (revision 962) @@ -23,7 +23,4 @@ char* outputfilename=NULL; char* lockname=NULL; - char* qmuinname=NULL; - char* qmuoutname=NULL; - char* qmuerrname=NULL; int numberofnodes; int waitonlock=0; @@ -62,7 +59,4 @@ outputfilename=argv[3]; lockname=argv[4]; - qmuinname=argv[5]; - qmuoutname=argv[6]; - qmuerrname=argv[7]; /*Open handle to data on disk: */ @@ -106,5 +100,5 @@ #ifdef _HAVE_DAKOTA_ - qmu(qmuinname,qmuoutname,qmuerrname,&fem,inputs,PrognosticAnalysisEnum(),NoneAnalysisEnum()); + Qmux(&fem,inputs,PrognosticAnalysisEnum(),NoneAnalysisEnum()); #else throw ErrorException(__FUNCT__," Dakota not present, cannot do qmu!"); Index: sm/trunk/src/c/parallel/qmu.cpp =================================================================== --- /issm/trunk/src/c/parallel/qmu.cpp (revision 961) +++ (revision ) @@ -1,126 +1,0 @@ -/*!\file: qmu.cpp - * \brief: wrapper to the Dakota capabilities. qmu fires up Dakota, and registers a Dakota Pluggin - * which will be in charge of running the solution sequences repeteadly, to garner statistics. - * - * This routine deals with running ISSM and Dakota in library mode. In library mode, Dakota does not - * run as an execuatble. Its capabilities are linked into the ISSM software. ISSM calls dakota routines - * directly from the dakota library. qmu.cpp is the code that is in charge of calling those routines. - * - * Dakota has its own way of running in parallel (for embarassingly parallel jobs). We do not want that, - * as ISSM knows exactly how to run "really parallel" jobs that use all CPUS. To bypass Dakota's parallelism, - * we overloaded the constructor for the parallel library (see the Dakota patch in the externalpackages/dakota - * directory). This overloaded constructor fires up Dakota serially on CPU 0 only! We take care of broadcasting - * to the other CPUS, hence ISSM is running in parallel, and Dakota serially on CPU0. - * - * Now, how does CPU 0 drive all other CPUS to carry out sensitivity analysese? By synchronizing its call to - * our ISSM cores (diagnostic_core, thermal_core, transient_core, etc ...) on CPU 0 with all other CPUS. - * This explains the structure of qmu.cpp, where cpu 0 runs Dakota, the Dakota pluggin fires up SpawnCore.cpp, - * while the other CPUS are waiting for a broadcast from CPU0, once they get it, they also fire up - * SpawnCore. In the end, SpawnCore is fired up on all CPUS, with CPU0 having Dakota inputs, that it will - * broacast to other CPUS. - * - * Now, how does dakota call the SpawnCore routine? The SpawnCore is embedded into the DakotaPlugin object - * which is derived from the Direct Interface Dakota objct. This is the only way to run Dakota in library - * mode (see their developper guide for more info). Dakota registers the DakotaPlugin object into its own - * database, and calls on the embedded SpawnCore from CPU0. - * - */ - -#ifdef HAVE_CONFIG_H -#include "config.h" -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif - -#undef __FUNCT__ -#define __FUNCT__ "qmu" - -#include "../objects/DakotaPlugin.h" -#include "../shared/shared.h" -#include "../include/macros.h" -#include "parallel.h" - - -#ifdef _HAVE_DAKOTA_ //only works if dakota library has been compiled in. - -#include "ParallelLibrary.H" -#include "ProblemDescDB.H" -#include "DakotaStrategy.H" -#include "DakotaModel.H" -#include "DakotaInterface.H" - - -void qmu(const char* dakota_input_file,const char* dakota_output_file,const char* dakota_error_file, FemModel* femmodels,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ - - - extern int my_rank; - int status=0; - Dakota::ModelLIter ml_iter; - - if(my_rank==0){ - - // Instantiate/initialize the parallel library and problem description - // database objects. - Dakota::ParallelLibrary parallel_lib("serial"); - Dakota::ProblemDescDB problem_db(parallel_lib); - - // Manage input file parsing, output redirection, and restart processing - // without a CommandLineHandler. This version relies on parsing of an - // input file. - problem_db.manage_inputs(dakota_input_file); - // specify_outputs_restart() is only necessary if specifying non-defaults - parallel_lib.specify_outputs_restart(dakota_output_file,dakota_error_file,NULL,NULL); - - // Instantiate the Strategy object (which instantiates all Model and - // Iterator objects) using the parsed information in problem_db. - Dakota::Strategy selected_strategy(problem_db); - - // convenience function for iterating over models and performing any - // interface plug-ins - Dakota::ModelList& models = problem_db.model_list(); - - for (ml_iter = models.begin(); ml_iter != models.end(); ml_iter++) { - - Dakota::Interface& interface = ml_iter->interface(); - - //set DB nodes to the existing Model specification - problem_db.set_db_model_nodes(ml_iter->model_id()); - - // Serial case: plug in derived Interface object without an analysisComm - interface.assign_rep(new SIM::DakotaPlugin(problem_db,femmodels,inputs,analysis_type,sub_analysis_type), false); - } - - // Execute the strategy - problem_db.lock(); // prevent run-time DB queries - selected_strategy.run_strategy(); - - //Warn other cpus that we are done running the dakota iterator - MPI_Bcast(&status,1,MPI_INT,0,MPI_COMM_WORLD); - - } - else{ - - for(;;){ - - SpawnCore(NULL,NULL,NULL,0,femmodels,inputs,analysis_type,sub_analysis_type,0); - - //Figure out if cpu 0 is done iterating - MPI_Bcast(&status,1,MPI_INT,0,MPI_COMM_WORLD); - - //yes! - if(status){ - break; //yes, we are done - } - else{ - /*We are not done: */ - - - } - } - - - } -} - -#endif //only works if dakota library has been compiled in. - Index: /issm/trunk/src/c/parallel/thermal.cpp =================================================================== --- /issm/trunk/src/c/parallel/thermal.cpp (revision 961) +++ /issm/trunk/src/c/parallel/thermal.cpp (revision 962) @@ -24,7 +24,4 @@ char* outputfilename=NULL; char* lockname=NULL; - char* qmuinname=NULL; - char* qmuoutname=NULL; - char* qmuerrname=NULL; int qmu_analysis=0; int numberofnodes; @@ -61,7 +58,4 @@ outputfilename=argv[3]; lockname=argv[4]; - qmuinname=argv[5]; - qmuoutname=argv[6]; - qmuerrname=argv[7]; /*Open handle to data on disk: */ @@ -114,5 +108,5 @@ #ifdef _HAVE_DAKOTA_ - qmu(qmuinname,qmuoutname,qmuerrname,&femmodels[0],inputs,ThermalAnalysisEnum(),NoneAnalysisEnum()); + Qmux(&femmodels[0],inputs,ThermalAnalysisEnum(),NoneAnalysisEnum()); #else throw ErrorException(__FUNCT__," Dakota not present, cannot do qmu!"); Index: /issm/trunk/src/c/parallel/transient.cpp =================================================================== --- /issm/trunk/src/c/parallel/transient.cpp (revision 961) +++ /issm/trunk/src/c/parallel/transient.cpp (revision 962) @@ -23,7 +23,4 @@ char* outputfilename=NULL; char* lockname=NULL; - char* qmuinname=NULL; - char* qmuoutname=NULL; - char* qmuerrname=NULL; int numberofnodes; int qmu_analysis=0; @@ -62,7 +59,4 @@ outputfilename=argv[3]; lockname=argv[4]; - qmuinname=argv[5]; - qmuoutname=argv[6]; - qmuerrname=argv[7]; /*Open handle to data on disk: */ @@ -124,5 +118,5 @@ #ifdef _HAVE_DAKOTA_ - qmu(qmuinname,qmuoutname,qmuerrname,&femmodels[0],inputs,TransientAnalysisEnum(),NoneAnalysisEnum()); + Qmux(&femmodels[0],inputs,TransientAnalysisEnum(),NoneAnalysisEnum()); #else throw ErrorException(__FUNCT__," Dakota not present, cannot do qmu!"); Index: /issm/trunk/src/mex/Makefile.am =================================================================== --- /issm/trunk/src/mex/Makefile.am (revision 961) +++ /issm/trunk/src/mex/Makefile.am (revision 962) @@ -29,4 +29,5 @@ PenaltySystemMatrices\ ProcessParams\ + Qmu\ Reduceloadfromgtof\ Reducematrixfromgtof\ @@ -48,5 +49,5 @@ endif -LDADD = $(MEXLIB) $(TRIANGLELIB) $(PETSCLIB) ../c/libISSM.a $(METISLIB) +LDADD = $(MEXLIB) $(TRIANGLELIB) $(PETSCLIB) ../c/libISSM.a $(DAKOTALIB) $(METISLIB) $(FLIBS) AM_LDFLAGS = $(MEXLINK) @@ -142,4 +143,7 @@ ProcessParams/ProcessParams.h +Qmu_SOURCES = Qmu/Qmu.cpp\ + Qmu/Qmu.h + PenaltySystemMatrices_SOURCES = PenaltySystemMatrices/PenaltySystemMatrices.cpp\ PenaltySystemMatrices/PenaltySystemMatrices.h Index: /issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp =================================================================== --- /issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp (revision 961) +++ /issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp (revision 962) @@ -33,5 +33,4 @@ /*Create elements, nodes and materials: */ CreateDataSets(&elements,&nodes,&materials,&constraints, &loads, ¶meters, model,MODEL); - /*Write output data: */ WriteData(ELEMENTS,elements,0,0,"DataSet",NULL); @@ -41,4 +40,5 @@ WriteData(MATERIALS,materials,0,0,"DataSet",NULL); WriteData(PARAMETERS,parameters,0,0,"DataSet",NULL); + /*Free ressources: */ Index: /issm/trunk/src/mex/Qmu/Qmu.cpp =================================================================== --- /issm/trunk/src/mex/Qmu/Qmu.cpp (revision 962) +++ /issm/trunk/src/mex/Qmu/Qmu.cpp (revision 962) @@ -0,0 +1,68 @@ +/*\file Qmu.c + *\brief: launch solution core using Dakota engine. + */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "./Qmu.h" + +void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){ + + /*diverse: */ + int noerr=1; + + /*input datasets: */ + mxArray* models=NULL; + mxArray* inputs=NULL; + int analysis_type; + int sub_analysis_type; + char* dakota_input_file=NULL; + char* dakota_output_file=NULL; + char* dakota_error_file=NULL; + + /* output datasets: none*/ + + #ifndef _HAVE_DAKOTA_ //only works if dakota library has been compiled in. + throw ErrorException(__FUNCT__," Dakota not available! Cannot carry out qmu analysis!"); + #endif + + /*Boot module: */ + MODULEBOOT(); + + /*checks on arguments on the matlab side: */ + CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&QmuUsage); + + /*Input datasets: */ + models=MODELS; + inputs=INPUTS; + + FetchData((void**)&analysis_type,NULL,NULL,mxGetField(PARAMETERS,0,"analysis_type"),"Integer",NULL); + FetchData((void**)&sub_analysis_type,NULL,NULL,mxGetField(PARAMETERS,0,"sub_analysis_type"),"Integer",NULL); + FetchData((void**)&dakota_input_file,NULL,NULL,mxGetField(PARAMETERS,0,"qmuinname"),"String",NULL); + FetchData((void**)&dakota_output_file,NULL,NULL,mxGetField(PARAMETERS,0,"qmuoutname"),"String",NULL); + FetchData((void**)&dakota_error_file,NULL,NULL,mxGetField(PARAMETERS,0,"qmuerrname"),"String",NULL); + + printf("%s %s %s\n",dakota_output_file,dakota_input_file,dakota_error_file); + + /*!Generate internal degree of freedom numbers: */ + Qmux(models,inputs,analysis_type,sub_analysis_type,dakota_input_file,dakota_output_file,dakota_error_file); + + /*end module: */ + MODULEEND(); + + /*Free ressources:*/ + xfree((void**)&dakota_input_file); + xfree((void**)&dakota_output_file); + xfree((void**)&dakota_error_file); + +} + +void QmuUsage(void){ + _printf_("\n"); + _printf_(" usage: %s(models,inputs,analysis_type,sub_analysis_type,dakota_input_file,dakota_outputa_file,dakota_error_file);\n",__FUNCT__); + _printf_("\n"); +} Index: /issm/trunk/src/mex/Qmu/Qmu.h =================================================================== --- /issm/trunk/src/mex/Qmu/Qmu.h (revision 962) +++ /issm/trunk/src/mex/Qmu/Qmu.h (revision 962) @@ -0,0 +1,33 @@ + +/* + Qmu.h +*/ + + +#ifndef _QMU_H +#define _QMU_H + +/* local prototypes: */ +void QmuUsage(void); + +#include "../../c/issm.h" + +#undef __FUNCT__ +#define __FUNCT__ "Qmu" + +/* serial input macros: */ +#define MODELS (mxArray*)prhs[0] +#define INPUTS (mxArray*)prhs[1] +#define PARAMETERS (mxArray*)prhs[2] + +/* serial output macros: */ + +/* serial arg counts: */ +#undef NLHS +#define NLHS 0 +#undef NRHS +#define NRHS 3 + + +#endif /* _QMU_H */ +