Changeset 16142

Timestamp:

09/16/13 22:05:45 (12 years ago)

Author:

Eric.Larour

Message:

CHG: added new type of matrix solver type, which allows us to make run time decisions
on the type of package we use to solve our matrix systems.

Location:

issm/trunk-jpl

Files:

• src/c/classes/FemModel.cpp (modified) (4 diffs)
• src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp (modified) (2 diffs)
• src/c/modules/ModelProcessorx/CreateDataSets.cpp (modified) (2 diffs)
• src/c/modules/ModelProcessorx/CreateParameters.cpp (modified) (2 diffs)
• src/c/modules/ModelProcessorx/ModelProcessorx.cpp (modified) (2 diffs)
• src/c/modules/ModelProcessorx/ModelProcessorx.h (modified) (1 diff)
• src/c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• src/c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• src/c/shared/Enum/StringToEnumx.cpp (modified) (1 diff)
• src/c/toolkits/issm/IssmMpiDenseMat.h (modified) (1 diff)
• src/c/toolkits/issm/IssmToolkitUtils.cpp (modified) (3 diffs)
• src/c/toolkits/issm/IssmToolkitUtils.h (modified) (1 diff)
• src/m/enum/EnumDefinitions.py (modified) (1 diff)
• src/m/enum/GslEnum.m (added)
• src/m/enum/MumpsEnum.m (added)
• src/m/solvers/issmsolver.m (modified) (1 diff)
• src/m/solvers/issmsolver.py (modified) (1 diff)
• test/NightlyRun/testad.m (modified) (1 diff)

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -189 +189 @@
         /*Open input file on cpu 0: */
         if(my_rank==0) IOMODEL = pfopen0(inputfilename ,"rb");
+        
+        /*Open toolkits file: */
+        toolkitsoptionsfid=pfopen(toolkitsfilename,"r");
 
         /*Initialize internal data: */
@@ -203 +206 @@
 
         /*create datasets for all analyses*/
-        ModelProcessorx(&this->elements,&this->nodes,&this->vertices,&this->materials,&this->constraints,&this->loads,&this->parameters,IOMODEL,rootpath,this->solution_type,nummodels,analyses);
+        ModelProcessorx(&this->elements,&this->nodes,&this->vertices,&this->materials,&this->constraints,&this->loads,&this->parameters,IOMODEL,toolkitsoptionsfid,rootpath,this->solution_type,nummodels,analyses);
 
         /*do the post-processing of the datasets to get an FemModel that can actually run analyses: */
@@ -227 +230 @@
         }
 
-        /*Close input file descriptors: */
+        /*Close input file and toolkits file descriptors: */
         if(my_rank==0) pfclose(IOMODEL,inputfilename);
+        pfclose(toolkitsoptionsfid,toolkitsfilename);
 
         /*Open output file once for all and add output file name and file descriptor to parameters*/
@@ -238 +242 @@
         this->parameters->AddObject(new StringParam(LockFileNameEnum,lockfilename));
 
-        /*Now, deal with toolkits options, which need to be put into the parameters dataset: */
-        toolkitsoptionsfid=pfopen(toolkitsfilename,"r");
-        ParseToolkitsOptionsx(this->parameters,toolkitsoptionsfid);
-        pfclose(toolkitsoptionsfid,toolkitsfilename);
-}
+        }
 /*}}}*/
 /*FUNCTION FemModel::OutputResults {{{*/

issm/trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp

@@ -21 +21 @@
         int*        indices=NULL;
         int         num_indices;
+        char* options=NULL;
 
         IssmDouble* xp=NULL;
@@ -116 +117 @@
                 /*initialize a placeholder to store solver pointers: {{{*/
                 GenericParam<Adolc_edf> *theAdolcEDF_p=new GenericParam<Adolc_edf>(AdolcParamEnum);
-                #ifdef _HAVE_GSL_
-                theAdolcEDF_p->GetParameterValue().myEDF_for_solverx_p=reg_ext_fct(EDF_for_solverx);
-                #endif
-                #ifdef _HAVE_MUMPS_
-                theAdolcEDF_p->GetParameterValue().myEDF_for_solverx_p=reg_ext_fct(mumpsSolveEDF);
-                #endif
+
+                /*Solver pointers depend on what type of solver we are implementing: */
+                options=OptionsFromAnalysis(parameters,DefaultAnalysisEnum); //options database is not filled in yet, use default.
+                ToolkitOptions::Init(options);
+
+                switch(IssmSolverTypeFromToolkitOptions()){
+                        case MumpsEnum:{
+                                #ifdef _HAVE_MUMPS_
+                                theAdolcEDF_p->GetParameterValue().myEDF_for_solverx_p=reg_ext_fct(mumpsSolveEDF);
+                                #else
+                                _error_("requesting mumps solver without MUMPS being compiled in!");
+                                #endif
+                                break;
+                                                        }
+                        case GslEnum: {
+                                #ifdef _HAVE_GSL_
+                                theAdolcEDF_p->GetParameterValue().myEDF_for_solverx_p=reg_ext_fct(EDF_for_solverx);
+                                #endif
+                            break;
+                                                  }
+                        default:
+                                _error_("solver type not supported yet!");
+                }
+
                 // to save some space:
                 // we know we won't use adolc inside of  the solver:

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateDataSets.cpp

@@ -14 +14 @@
 #include "./ModelProcessorx.h"
 
-void CreateDataSets(Elements** pelements,Nodes** pnodes, Vertices** pvertices, Materials** pmaterials, Constraints** pconstraints, Loads** ploads,Parameters** pparameters,IoModel* iomodel,char* rootpath,const int solution_type,const int analysis_type,const int nummodels,int analysis_counter){
+void CreateDataSets(Elements** pelements,Nodes** pnodes, Vertices** pvertices, Materials** pmaterials, Constraints** pconstraints, Loads** ploads,Parameters** pparameters,IoModel* iomodel,FILE* toolkitfile,char* rootpath,const int solution_type,const int analysis_type,const int nummodels,int analysis_counter){
 
         Elements   *elements   = NULL;
@@ -178 +178 @@
 
         /*Generate objects that are not dependent on any analysis_type: */
-        CreateParameters(pparameters,iomodel,rootpath,solution_type,analysis_type,analysis_counter);
+        CreateParameters(pparameters,iomodel,rootpath,toolkitfile,solution_type,analysis_type,analysis_counter);
 
         /*Update Elements in case we are running a transient solution: */

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

@@ -13 +13 @@
 #include "../../shared/shared.h"
 #include "../MeshPartitionx/MeshPartitionx.h"
+#include "../ParseToolkitsOptionsx/ParseToolkitsOptionsx.h"
 #include "./ModelProcessorx.h"
 
-void CreateParameters(Parameters** pparameters,IoModel* iomodel,char* rootpath,const int solution_type,int analysis_type,int analysis_counter){
+void CreateParameters(Parameters** pparameters,IoModel* iomodel,char* rootpath,FILE* toolkitsoptionsfid,const int solution_type,int analysis_type,int analysis_counter){
 
         int         i,j,m,k;
@@ -244 +245 @@
         CreateParametersDakota(&parameters,iomodel,rootpath,solution_type,analysis_type);
         #endif
+        
+        /*Now, deal with toolkits options, which need to be put into the parameters dataset: */
+        ParseToolkitsOptionsx(parameters,toolkitsoptionsfid);
 
         #ifdef _HAVE_ADOLC_

issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp

@@ -13 +13 @@
 #include "./ModelProcessorx.h"
 
-void ModelProcessorx(Elements** pelements, Nodes** pnodes, Vertices** pvertices, Materials** pmaterials, Constraints** pconstraints, Loads** ploads, Parameters** pparameters, FILE* IOMODEL,char* rootpath,const int solution_type,const int nummodels,const int* analysis_type_list){
+void ModelProcessorx(Elements** pelements, Nodes** pnodes, Vertices** pvertices, Materials** pmaterials, Constraints** pconstraints, Loads** ploads, Parameters** pparameters, FILE* IOMODEL,FILE* toolkitfile, char* rootpath,const int solution_type,const int nummodels,const int* analysis_type_list){
 
         int   i,analysis_type,verbose;
@@ -65 +65 @@
 
                 if(VerboseMProcessor()) _printf0_("   creating datasets for analysis " << EnumToStringx(analysis_type) << "\n");
-                CreateDataSets(&elements,&nodes,&vertices,&materials,&constraints,&loads,&parameters,iomodel,rootpath,solution_type,analysis_type,nummodels,i);
+                CreateDataSets(&elements,&nodes,&vertices,&materials,&constraints,&loads,&parameters,iomodel,toolkitfile,rootpath,solution_type,analysis_type,nummodels,i);
         }
         if(VerboseMProcessor()) _printf0_("   done with model processor \n");

issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h

@@ -10 +10 @@
 #include "../../classes/classes.h"
 
-void ModelProcessorx(Elements** pelements, Nodes** pnodes, Vertices** pvertices, Materials** pmaterials, Constraints** pconstraints, Loads** ploads, Parameters** pparameters, FILE* iomodel_handle,char* rootpath,const int solution_type,const int nummodels,const int* analysis_type_listh);
+void ModelProcessorx(Elements** pelements, Nodes** pnodes, Vertices** pvertices, Materials** pmaterials, Constraints** pconstraints, Loads** ploads, Parameters** pparameters, FILE* iomodel_handle,FILE* toolkitfile, char* rootpath,const int solution_type,const int nummodels,const int* analysis_type_listh);
 
 /*Creation of fem datasets: general drivers*/
-void CreateDataSets(Elements** pelements,Nodes** pnodes,Vertices** pvertices, Materials** pmaterials, Constraints** pconstraints, Loads** ploads,Parameters** pparameters,IoModel* iomodel,char* rootpath,const int solution_type,int analysis_type,const int nummodels,int analysis_counter);
+void CreateDataSets(Elements** pelements,Nodes** pnodes,Vertices** pvertices, Materials** pmaterials, Constraints** pconstraints, Loads** ploads,Parameters** pparameters,IoModel* iomodel,FILE* toolkitfile,char* rootpath,const int solution_type,int analysis_type,const int nummodels,int analysis_counter);
 void CreateElementsVerticesAndMaterials(Elements** pelements,Vertices** pvertices,Materials** pmaterials, IoModel* iomodel,const int nummodels);
-void CreateParameters(Parameters** pparameters,IoModel* iomodel,char* rootpath,const int solution_type,int analysis_type,int analysis_counter);
+void CreateParameters(Parameters** pparameters,IoModel* iomodel,char* rootpath,FILE* toolkitfile,const int solution_type,int analysis_type,int analysis_counter);
 void CreateParametersAutodiff(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type);
 void CreateParametersControl(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type);

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -608 +608 @@
         SeqEnum,
         MpiEnum,
+        MumpsEnum,
+        GslEnum,
         /*}}}*/
         /*Options{{{*/

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -582 +582 @@
                 case SeqEnum : return "Seq";
                 case MpiEnum : return "Mpi";
+                case MumpsEnum : return "Mumps";
+                case GslEnum : return "Gsl";
                 case OptionEnum : return "Option";
                 case GenericOptionEnum : return "GenericOption";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -594 +594 @@
               else if (strcmp(name,"Seq")==0) return SeqEnum;
               else if (strcmp(name,"Mpi")==0) return MpiEnum;
+              else if (strcmp(name,"Mumps")==0) return MumpsEnum;
+              else if (strcmp(name,"Gsl")==0) return GslEnum;
               else if (strcmp(name,"Option")==0) return OptionEnum;
               else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum;

issm/trunk-jpl/src/c/toolkits/issm/IssmMpiDenseMat.h

@@ -513 +513 @@
                         pf=(IssmMpiVec<IssmDouble>*)pfin;
 
-                        /*Initialize output: */
-                        uf=pf->Duplicate();
-
-                        /*Let's try and use the MUMPS solver here: */
-                        #ifdef _HAVE_MUMPS_
-                        MpiDenseMumpsSolve(/*output*/ uf->vector,uf->M,uf->m, /*stiffness matrix:*/ this->matrix,this->M,this->N,this->m, /*right hand side load vector: */ pf->vector,pf->M,pf->m,parameters);
-                        #else
-                        _error_("IssmMpiDenseMat solver requires MUMPS solver");
-                        #endif
-                        return (IssmAbsVec<IssmDouble>*)uf;
+                                                
+                        switch(IssmSolverTypeFromToolkitOptions()){
+                                case MumpsEnum: {
+                                        /*Initialize output: */
+                                        uf=pf->Duplicate();
+                                        #ifdef _HAVE_MUMPS_
+                                        MpiDenseMumpsSolve(/*output*/ uf->vector,uf->M,uf->m, /*stiffness matrix:*/ this->matrix,this->M,this->N,this->m, /*right hand side load vector: */ pf->vector,pf->M,pf->m,parameters);
+                                        #else
+                                        _error_("IssmMpiDenseMat solver requires MUMPS solver");
+                                        #endif
+                                        return (IssmAbsVec<IssmDouble>*)uf;
+                                        break;
+                                                                }
+                                case GslEnum: {
+
+                                        IssmDouble* x=NULL;
+                                        #ifdef _HAVE_GSL_
+
+                                        DenseGslSolve(/*output*/ &x,/*stiffness matrix:*/ this->matrix,this->M,this->N, /*right hand side load vector: */ pf->vector,pf->M,parameters);
+
+                                        uf=new IssmMpiVec<IssmDouble>(x,this->N); xDelete(x);
+
+                                        return (IssmAbsVec<IssmDouble>*)uf;
+                                        break;
+                                        #else
+                                        _error_("GSL support not compiled in!");
+                                        #endif
+                                                          }
+                                default:
+                                        _error_("solver type not supported yet!");
+                        }
 
                 }/*}}}*/

issm/trunk-jpl/src/c/toolkits/issm/IssmToolkitUtils.cpp

@@ -36 +36 @@
         mat_type=ToolkitOptions::GetToolkitOptionValue("mat_type");
 
-        if ((strcmp(mat_type,"mpidense")==0) || (strcmp(mat_type,"dense")==0)){
-                if (isparallel) mat_type_enum=MpiDenseEnum;
+        if (strcmp(mat_type,"mpidense")==0){
+                mat_type_enum=MpiDenseEnum;
+        }
+        else if (strcmp(mat_type,"dense")==0){
+                if (isparallel) _error_("Dense matrix type not supported for parallel runs with num_procs>1");
                 else mat_type_enum=DenseEnum;
         }
@@ -62 +65 @@
          *we try and stick with the Petsc vector types: */
         vec_type=ToolkitOptions::GetToolkitOptionValue("vec_type");
-
-        if ((strcmp(vec_type,"mpi")==0) || (strcmp(vec_type,"seq")==0)){
-                if (isparallel) vec_type_enum=MpiEnum;
+        
+        if (strcmp(vec_type,"mpi")==0){
+                vec_type_enum=MpiEnum;
+        }
+        else if (strcmp(vec_type,"seq")==0){
+                if (isparallel) _error_("Dense vector type not supported for parallel runs with num_procs>1");
                 else vec_type_enum=SeqEnum;
         }
@@ -75 +81 @@
         return vec_type_enum;
 } /*}}}*/  
+int IssmSolverTypeFromToolkitOptions(void){ /*{{{*/
+
+        char* solver_type=NULL;
+        int   solver_type_enum;
+        int   num_procs=0;
+        bool  isparallel=false;
+
+        /*first, figure out if we are running in parallel: */
+        num_procs=IssmComm::GetSize();
+        if(num_procs>1)isparallel=true;
+
+        /*retrieve solver type as a string, from the Toolkits Options database, similar to what Petsc does. Actually, 
+         *we try and stick with the Petsc vector types: */
+        solver_type=ToolkitOptions::GetToolkitOptionValue("solver_type");
+
+        if (strcmp(solver_type,"mumps")==0){
+                solver_type_enum=MumpsEnum;
+        }
+        else if (strcmp(solver_type,"gsl")==0){
+                if (isparallel) _error_("Gsl solver type not supported for parallel runs with num_procs>1");
+                else solver_type_enum=GslEnum;
+        }
+        else _error_("solver type not supported yet!");
+
+        /*free ressources: */
+        xDelete<char>(solver_type);
+
+        /*return: */
+        return solver_type_enum;
+} /*}}}*/  

issm/trunk-jpl/src/c/toolkits/issm/IssmToolkitUtils.h

@@ -18 +18 @@
 int IssmMatTypeFromToolkitOptions(void);
 int IssmVecTypeFromToolkitOptions(void);
+int IssmSolverTypeFromToolkitOptions(void);
 
 #endif //#ifndef _ISSMMAT_H_

issm/trunk-jpl/src/m/enum/EnumDefinitions.py

@@ -574 +574 @@
 def SeqEnum(): return StringToEnum("Seq")[0]
 def MpiEnum(): return StringToEnum("Mpi")[0]
+def MumpsEnum(): return StringToEnum("Mumps")[0]
+def GslEnum(): return StringToEnum("Gsl")[0]
 def OptionEnum(): return StringToEnum("Option")[0]
 def GenericOptionEnum(): return StringToEnum("GenericOption")[0]

issm/trunk-jpl/src/m/solvers/issmsolver.m

@@ -13 +13 @@
 issmoptions.mat_type=getfieldvalue(options,'mat_type','mpidense');
 issmoptions.vec_type=getfieldvalue(options,'vec_type','mpi');
+issmoptions.solver_type=getfieldvalue(options,'solver_type','mumps');

issm/trunk-jpl/src/m/solvers/issmsolver.py

@@ -9 +9 @@
         arguments=pairoptions(*args) 
         
-        options=[['toolkit','issm'],['mat_type','mpidense'],['vec_type','mpi']];
+        options=[['toolkit','issm'],['mat_type','mpidense'],['vec_type','mpi'],['solver_type','mumps']];
 
         #now, go through our arguments, and write over default options.

issm/trunk-jpl/test/NightlyRun/testad.m

@@ -1 +1 @@
 %test reverse scalar vs forward vectorial drivers in ADOLC, using the test3009 setup, equivalent to test109 setup.
-md=triangle(model(),'../Exp/Square.exp',100000.);
+md=triangle(model(),'../Exp/Square.exp',500000.);
 md=setmask(md,'all','');
 md=parameterize(md,'../Par/SquareShelfConstrained.par');
 md=setflowequation(md,'SSA','all');
-md.cluster=generic('name',oshostname(),'np',1);
+md.cluster=generic('name',oshostname(),'np',2);
 %md.debug.valgrind=true;
 %md.verbose=verbose('11111111')
 md.toolkits.DefaultAnalysis=issmsolver();
+%md.toolkits.DefaultAnalysis.solver_type='gsl';
 
 md.autodiff.isautodiff=true;