Changeset 659

Timestamp:

05/31/09 00:20:31 (16 years ago)

Author:

Eric.Larour

Message:

Finished Dakota responses in parallel

Location:

issm/trunk/src

Files:

• c/DataSet/DataSet.cpp (modified) (1 diff)
• c/DataSet/DataSet.h (modified) (1 diff)
• c/parallel/DakotaResponses.cpp (modified) (2 diffs)
• c/parallel/ProcessResults.cpp (modified) (9 diffs)
• c/parallel/SpawnCore.cpp (modified) (4 diffs)
• c/parallel/diagnostic.cpp (modified) (5 diffs)
• c/parallel/parallel.h (modified) (1 diff)
• c/parallel/qmu.cpp (modified) (3 diffs)
• c/parallel/thermal.cpp (modified) (3 diffs)
• m/classes/public/BuildQueueingScriptGeneric.m (modified) (1 diff)
• m/classes/public/LaunchQueueJobGeneric.m (modified) (2 diffs)
• m/classes/public/marshall.m (modified) (1 diff)
• m/solutions/dakota/dakota_in_write.m (modified) (1 diff)
• m/solutions/dakota/qmuin.m (modified) (2 diffs)

issm/trunk/src/c/DataSet/DataSet.cpp

@@ -317 +317 @@
 }
 
+int   DataSet::FindResult(void* pvalue, char* name){
+
+        /*Go through a dataset, and find a Result* object 
+         *which field name is "name" : */
+        
+        vector<Object*>::iterator object;
+        Result* result=NULL;
+
+        int found=0;
+
+        for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+                /*Find param type objects: */
+                if((*object)->Enum()==ResultEnum()){
+
+                        /*Ok, this object is a result, recover it and ask which name it has: */
+                        result=(Result*)(*object);
+
+                        if (strcmp(result->GetFieldName(),name)==0){
+                                /*Ok, this is the one! Recover the value of this result: */
+                                double** field=NULL;
+                                field=(double**)pvalue; //not very safe, but hey!
+                                result->GetField(field);
+                                found=1;
+                                break;
+                        }
+                }
+        }
+        return found;
+}
+
+
 Object*   DataSet::FindParamObject(char* name){
 

issm/trunk/src/c/DataSet/DataSet.h

@@ -81 +81 @@
                 int   DeleteObject(Object* object);
                 void  ComputePressure(Vec p_g);
+                int   FindResult(void* pvalue, char* name);
 
 };

issm/trunk/src/c/parallel/DakotaResponses.cpp

@@ -14 +14 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "DakotaResponses"
-void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,DataSet* results,int analysis_type,int sub_analysis_type){
+void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,FemModel* femmodels,DataSet* results,int analysis_type,int sub_analysis_type){
 
-        int i;
+        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;i<numresponses;i++){
@@ -25 +37 @@
                 //'min_vx' 'max_vx' 'max_abs_vx' 'min_vy' 'max_vy' 'max_abs_vy' 'min_vel' 'max_vel'
 
-                if(strcmp(response_descriptor,"min_vel")){
+                if(strcmp(response_descriptor,"min_vel")==0){
+                        double* vel=NULL;
+                        double min_vel=0;
+                
+                        found=results->FindResult((void*)&vel,"vel");
+                        if(!found)throw ErrorException(__FUNCT__," could not find vel to compute min_vel");
 
+                        min_vel=vel[0];
+                        for(j=1;j<numberofnodes;j++){
+                                if (vel[j]<min_vel)min_vel=vel[j];
+                        }
 
+                        if(my_rank==0)responses[i]=min_vel;
+                        
+                }
+                else if(strcmp(response_descriptor,"max_vel")==0){
+                        double* vel=NULL;
+                        double max_vel=0;
+
+                        found=results->FindResult((void*)&vel,"vel");
+                        if(!found)throw ErrorException(__FUNCT__," could not find vel to compute max_vel");
+
+                        max_vel=vel[0];
+                        for(j=1;j<numberofnodes;j++){
+                                if (vel[j]>max_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;j<numberofnodes;j++){
+                                if (vx[j]<min_vx)min_vx=vx[j];
+                        }
+                        if(my_rank==0)responses[i]=min_vx;
+                }
+                else if(strcmp(response_descriptor,"max_vx")==0){
+                        double* vx=NULL;
+                        double max_vx=0;
+                        
+                        found=results->FindResult((void*)&vx,"vx");
+                        if(!found)throw ErrorException(__FUNCT__," could not find vx to compute max_vx");
+
+                        max_vx=vx[0];
+                        for(j=1;j<numberofnodes;j++){
+                                if (vx[j]>max_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;j<numberofnodes;j++){
+                                if (fabs(vx[j])>max_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;j<numberofnodes;j++){
+                                if (vy[j]<min_vy)min_vy=vy[j];
+                        }
+                        if(my_rank==0)responses[i]=min_vy;
+                }
+                else if(strcmp(response_descriptor,"max_vy")==0){
+                        double* vy=NULL;
+                        double max_vy=0;
+                        
+                        found=results->FindResult((void*)&vy,"vy");
+                        if(!found)throw ErrorException(__FUNCT__," could not find vy to compute max_vy");
+
+                        max_vy=vy[0];
+                        for(j=1;j<numberofnodes;j++){
+                                if (vy[j]>max_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;j<numberofnodes;j++){
+                                if (fabs(vy[j])>max_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;j<numberofnodes;j++){
+                                if (vz[j]<min_vz)min_vz=vz[j];
+                        }
+                        if(my_rank==0)responses[i]=min_vz;
+                }
+                else if(strcmp(response_descriptor,"max_vz")==0){
+                        double* vz=NULL;
+                        double max_vz=0;
+                        
+                        found=results->FindResult((void*)&vz,"vz");
+                        if(!found)throw ErrorException(__FUNCT__," could not find vz to compute max_vz");
+
+                        max_vz=vz[0];
+                        for(j=1;j<numberofnodes;j++){
+                                if (vz[j]>max_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;j<numberofnodes;j++){
+                                if (fabs(vz[j])>max_abs_vz)max_abs_vz=fabs(vz[j]);
+                        }
+                        if(my_rank==0)responses[i]=max_abs_vz;
                 }
                 else throw ErrorException(__FUNCT__,exprintf("%s%s%s"," response descriptor : ",response_descriptor," not supported yet!"));

issm/trunk/src/c/parallel/ProcessResults.cpp

@@ -21 +21 @@
 #include "../shared/shared.h"
 
-void ProcessResults(DataSet** pnewresults,DataSet* results,FemModel* fems,int analysis_type){
+void ProcessResults(DataSet** presults,FemModel* fems,int analysis_type){
 
         int i,n;
         Result* result=NULL;
         Result* newresult=NULL;
+        
+        /*input: */
+        DataSet* results=NULL;
 
         /*output: */
@@ -47 +50 @@
         double* vy=NULL;
         double* vz=NULL;
+        double* vel=NULL;
         Vec     p_g=NULL;
         double* p_g_serial=NULL;
@@ -55 +59 @@
         int numberofnodes;
 
+        /*recover input results: */
+        results=*presults;
+
         /*Initialize new results: */
         newresults=new DataSet(ResultsEnum());
                 
-
         /*Recover femmodels first: */
         if(analysis_type==DiagnosticAnalysisEnum()){
@@ -101 +107 @@
                                 vy=(double*)xmalloc(numberofnodes*sizeof(double));
                                 vz=(double*)xmalloc(numberofnodes*sizeof(double)); 
+                                vel=(double*)xmalloc(numberofnodes*sizeof(double)); 
 
                                 for(i=0;i<numberofnodes;i++){
@@ -106 +113 @@
                                         vy[i]=u_g_serial[2*(int)partition[i]+1]*yts;
                                         vz[i]=0;
+                                        vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
                                 }
+
                         }
                         else{
@@ -116 +125 @@
                                 vy=(double*)xmalloc(numberofnodes*sizeof(double));
                                 vz=(double*)xmalloc(numberofnodes*sizeof(double));
+                                vel=(double*)xmalloc(numberofnodes*sizeof(double));
                                 for(i=0;i<numberofnodes;i++){
                                         vx[i]=u_g_serial[4*(int)partition[i]+0]*yts;
                                         vy[i]=u_g_serial[4*(int)partition[i]+1]*yts;
                                         vz[i]=u_g_serial[4*(int)partition[i]+2]*yts;
+                                        vel[i]=sqrt(pow(vx[i],2)+pow(vy[i],2)+pow(vz[i],2));
                                 }
                         }
@@ -133 +144 @@
                         newresults->AddObject(newresult);
 
+                        newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vel",vel,numberofnodes);
+                        newresults->AddObject(newresult);
+
                         /*do some cleanup: */
                         xfree((void**)&u_g_serial);
@@ -138 +152 @@
                 }
                 else if(strcmp(result->GetFieldName(),"p_g")==0){
-
                         /*easy, p_g is of size numberofnodes, on 1 dof, just repartition: */
                         result->GetField(&p_g);
@@ -174 +187 @@
         }
 
+        /*Delete results: */
+        delete results;
+
+
         /*Assign output pointers:*/
-        *pnewresults=newresults;
-
+        *presults=newresults;
 }
 

issm/trunk/src/c/parallel/SpawnCore.cpp

@@ -24 +24 @@
         
         /*output from core solutions: */
-        Vec u_g=NULL;
-        Vec p_g=NULL;
+        DataSet* results=NULL;
 
         char** responses_descriptors=NULL;
@@ -34 +33 @@
         double* qmu_part=NULL;
         int     qmu_npart;
-        DataSet* results=NULL;
 
         extern int my_rank;
@@ -91 +89 @@
         #endif
 
+        _printf_("initialize results:\n");
+        results=new DataSet(ResultsEnum());
+
         /*Modify core inputs to reflect the dakota variables inputs: */
         //inputs->UpdateFromDakota(variables,variables_descriptors,numvariables,qmu_part,qmu_npart);
@@ -105 +106 @@
         }
 
+        /*Now process the outputs, before computing the dakota responses: */
+        _printf_("process results:\n");
+        ProcessResults(&results,femmodels,analysis_type); 
+
         /*compute responses on cpu 0: dummy for now! */
-        DakotaResponses(responses,responses_descriptors,numresponses,results,analysis_type,sub_analysis_type);
+        _printf_("compute dakota responses:\n");
+        DakotaResponses(responses,responses_descriptors,numresponses,femmodels,results,analysis_type,sub_analysis_type);
 
         /*Free ressources:*/
-        
-        //vectors
-        VecFree(&u_g);
-        VecFree(&p_g);
-        
+        delete results;
+
         //variables only on cpu != 0
         if(my_rank!=0){

issm/trunk/src/c/parallel/diagnostic.cpp

@@ -23 +23 @@
         char* outputfilename=NULL;
         char* lockname=NULL;
-        char* qmuname=NULL;
+        char* qmuinname=NULL;
+        char* qmuoutname=NULL;
+        char* qmuerrname=NULL;
         int   numberofnodes;
         int   qmu_analysis=0;
@@ -32 +34 @@
         /*Results: */
         DataSet* results=NULL;
-        DataSet* newresults=NULL;
         
         ParameterInputs* inputs=NULL;
@@ -56 +57 @@
         outputfilename=argv[3];
         lockname=argv[4];
-        qmuname=argv[5];
+        qmuinname=argv[5];
+        qmuoutname=argv[6];
+        qmuerrname=argv[7];
 
         /*Open handle to data on disk: */
@@ -101 +104 @@
                 _printf_("calling qmu analysis on diagnostic core:\n");
                 
-                //qmu(qmuname,&femmodels[0],inputs,DiagnosticAnalysisEnum(),NoneAnalysisEnum());
+                qmu(qmuinname,qmuoutname,qmuerrname,&femmodels[0],inputs,DiagnosticAnalysisEnum(),NoneAnalysisEnum());
         }
 
         _printf_("process results:\n");
-        ProcessResults(&newresults,results,&femmodels[0],DiagnosticAnalysisEnum()); delete results;
+        ProcessResults(&results,&femmodels[0],DiagnosticAnalysisEnum());
         
         _printf_("write results to disk:\n");
-        OutputResults(newresults,outputfilename);
+        OutputResults(results,outputfilename);
 
         _printf_("write lock file:\n");
@@ -117 +120 @@
                 
         _printf_("closing MPI and Petsc\n");
-        MPI_Barrier(MPI_COMM_WORLD);
-
-        /*Close MPI libraries: */
         PetscFinalize(); 
+        
 
         /*end module: */

issm/trunk/src/c/parallel/parallel.h

@@ -36 +36 @@
 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,FemModel* femmodels,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
+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);
-void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,DataSet* results,int analysis_type,int sub_analysis_type);
-void ProcessResults(DataSet** pnewresults,DataSet* results,FemModel* fems,int analysis_type);
+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);
 
 #endif

issm/trunk/src/c/parallel/qmu.cpp

@@ -23 +23 @@
 #include "parallel.h"
 
-void qmu(const char* dakota_input_file,FemModel* femmodels,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
+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;
@@ -30 +31 @@
 
         if(my_rank==0){
-
+        
                 // Instantiate/initialize the parallel library and problem description
                 // database objects.
@@ -41 +42 @@
                 problem_db.manage_inputs(dakota_input_file);
                 // specify_outputs_restart() is only necessary if specifying non-defaults
-                //parallel_lib.specify_outputs_restart(NULL, NULL, NULL, NULL);
+                parallel_lib.specify_outputs_restart(dakota_output_file,dakota_error_file,NULL, NULL);
 
                 // Instantiate the Strategy object (which instantiates all Model and

issm/trunk/src/c/parallel/thermal.cpp

@@ -35 +35 @@
         /*Results: */
         DataSet* results=NULL;
-        DataSet* newresults=NULL;
         
         ParameterInputs* inputs=NULL;
@@ -89 +88 @@
         param=(Param*)femmodels[1].parameters->FindParamObject("p_g");
         femmodels[1].parameters->DeleteObject((Object*)param);
+
         
         _printf_("call computational core:\n");
@@ -94 +94 @@
 
         _printf_("process results:\n");
-        ProcessResults(&newresults,results,&femmodels[0],ThermalAnalysisEnum()); delete results;
+        ProcessResults(&results,&femmodels[0],ThermalAnalysisEnum());
 
         _printf_("write results to disk:\n");
-        OutputResults(newresults,outputfilename);
+        OutputResults(results,outputfilename);
 
         _printf_("write lock file:\n");

issm/trunk/src/m/classes/public/BuildQueueingScriptGeneric.m

@@ -29 +29 @@
 end
 
-fprintf(fid,' %s %s.bin %s.outbin %s.lock qmu.in 2> %s.errlog >%s.outlog & ',executionpath,md.name,md.name,md.name,md.name,md.name);
+fprintf(fid,' %s %s.bin %s.outbin %s.lock %s.qmu.in %s.qmu.out %s.qmu.err 2> %s.errlog >%s.outlog & ',executionpath,md.name,md.name,md.name,md.name,md.name,md.name,md.name,md.name);
 fclose(fid);

issm/trunk/src/m/classes/public/LaunchQueueJobGeneric.m

@@ -21 +21 @@
 if strcmpi(hostname,md.cluster),
         if md.qmu_analysis,
-                system(['cp ' md.name '.bin' ' ' md.name '.queue qmu/qmu.in ' ' ' executionpath]);
+                system(['cp ' md.name '.bin' ' ' md.name '.queue qmu/' md.name '.qmu.in ' ' ' executionpath]);
         else
                 system(['cp ' md.name '.bin' ' ' md.name '.queue' ' ' executionpath]);
@@ -27 +27 @@
 else
         if md.qmu_analysis,
-                system(['scp ' md.name '.bin' ' ' md.name '.queue qmu/qmu.in ' ' ' md.cluster ':' executionpath]);
+                system(['scp ' md.name '.bin' ' ' md.name '.queue qmu/' md.name '.qmu.in ' ' ' md.cluster ':' executionpath]);
         else
                 system(['scp ' md.name '.bin' ' ' md.name '.queue' ' ' md.cluster ':' executionpath]);

issm/trunk/src/m/classes/public/marshall.m

@@ -164 +164 @@
                 WriteData(fid,responses(i).descriptor,'String',['descriptor' num2str(i-1)]);
         end
+        WriteData(fid,md.npart,'Integer','npart');
 end
         

issm/trunk/src/m/solutions/dakota/dakota_in_write.m

@@ -26 +26 @@
 [pathstr,name,ext,versn] = fileparts(filei);
 if isempty(ext)
-    ext='.in';
+    ext='.qmu.in';
 end
 filei2=fullfile(pathstr,[name ext versn]);

issm/trunk/src/m/solutions/dakota/qmuin.m

@@ -8 +8 @@
 qmudir ='qmu';
 %  qmufile can not be changed unless cielo_ice_script.sh is also changed
-qmufile='qmu';
+qmufile=md.name;
 ivar   =1;
 iresp  =1;
@@ -75 +75 @@
     md.qmu_params(iparams).analysis_driver=[ISSM_DIR '/src/m/solutions/dakota/cielo_ice_script.sh'];
 end
+
 dakota_in_data(md.qmu_method(imethod),md.variables(ivar),md.responses(iresp),md.qmu_params(iparams),qmufile,package,md);
 
-system(['rm -rf qmu.m']);
+
+system(['rm -rf ' md.name '.m']);
 cd ../