Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 16141) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 16142) @@ -189,4 +189,7 @@ /*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,5 +206,5 @@ /*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,6 +230,7 @@ } - /*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,9 +242,5 @@ 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 {{{*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp (revision 16141) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp (revision 16142) @@ -21,4 +21,5 @@ int* indices=NULL; int num_indices; + char* options=NULL; IssmDouble* xp=NULL; @@ -116,10 +117,28 @@ /*initialize a placeholder to store solver pointers: {{{*/ GenericParam *theAdolcEDF_p=new GenericParam(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: Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 16141) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 16142) @@ -14,5 +14,5 @@ #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,5 +178,5 @@ /*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: */ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 16141) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 16142) @@ -13,7 +13,8 @@ #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,4 +245,7 @@ CreateParametersDakota(¶meters,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_ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 16141) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 16142) @@ -13,5 +13,5 @@ #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,5 +65,5 @@ if(VerboseMProcessor()) _printf0_(" creating datasets for analysis " << EnumToStringx(analysis_type) << "\n"); - CreateDataSets(&elements,&nodes,&vertices,&materials,&constraints,&loads,¶meters,iomodel,rootpath,solution_type,analysis_type,nummodels,i); + CreateDataSets(&elements,&nodes,&vertices,&materials,&constraints,&loads,¶meters,iomodel,toolkitfile,rootpath,solution_type,analysis_type,nummodels,i); } if(VerboseMProcessor()) _printf0_(" done with model processor \n"); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 16141) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 16142) @@ -10,10 +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); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16141) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16142) @@ -608,4 +608,6 @@ SeqEnum, MpiEnum, + MumpsEnum, + GslEnum, /*}}}*/ /*Options{{{*/ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16141) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16142) @@ -582,4 +582,6 @@ case SeqEnum : return "Seq"; case MpiEnum : return "Mpi"; + case MumpsEnum : return "Mumps"; + case GslEnum : return "Gsl"; case OptionEnum : return "Option"; case GenericOptionEnum : return "GenericOption"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16141) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16142) @@ -594,4 +594,6 @@ 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; Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiDenseMat.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiDenseMat.h (revision 16141) +++ /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiDenseMat.h (revision 16142) @@ -513,14 +513,35 @@ pf=(IssmMpiVec*)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*)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*)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(x,this->N); xDelete(x); + + return (IssmAbsVec*)uf; + break; + #else + _error_("GSL support not compiled in!"); + #endif + } + default: + _error_("solver type not supported yet!"); + } }/*}}}*/ Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmToolkitUtils.cpp =================================================================== --- /issm/trunk-jpl/src/c/toolkits/issm/IssmToolkitUtils.cpp (revision 16141) +++ /issm/trunk-jpl/src/c/toolkits/issm/IssmToolkitUtils.cpp (revision 16142) @@ -36,6 +36,9 @@ 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,7 +65,10 @@ *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,2 +81,32 @@ 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(solver_type); + + /*return: */ + return solver_type_enum; +} /*}}}*/ Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmToolkitUtils.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/issm/IssmToolkitUtils.h (revision 16141) +++ /issm/trunk-jpl/src/c/toolkits/issm/IssmToolkitUtils.h (revision 16142) @@ -18,4 +18,5 @@ int IssmMatTypeFromToolkitOptions(void); int IssmVecTypeFromToolkitOptions(void); +int IssmSolverTypeFromToolkitOptions(void); #endif //#ifndef _ISSMMAT_H_ Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16141) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16142) @@ -574,4 +574,6 @@ 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] Index: /issm/trunk-jpl/src/m/enum/GslEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/GslEnum.m (revision 16142) +++ /issm/trunk-jpl/src/m/enum/GslEnum.m (revision 16142) @@ -0,0 +1,11 @@ +function macro=GslEnum() +%GSLENUM - Enum of Gsl +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=GslEnum() + +macro=StringToEnum('Gsl'); Index: /issm/trunk-jpl/src/m/enum/MumpsEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/MumpsEnum.m (revision 16142) +++ /issm/trunk-jpl/src/m/enum/MumpsEnum.m (revision 16142) @@ -0,0 +1,11 @@ +function macro=MumpsEnum() +%MUMPSENUM - Enum of Mumps +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=MumpsEnum() + +macro=StringToEnum('Mumps'); Index: /issm/trunk-jpl/src/m/solvers/issmsolver.m =================================================================== --- /issm/trunk-jpl/src/m/solvers/issmsolver.m (revision 16141) +++ /issm/trunk-jpl/src/m/solvers/issmsolver.m (revision 16142) @@ -13,2 +13,3 @@ issmoptions.mat_type=getfieldvalue(options,'mat_type','mpidense'); issmoptions.vec_type=getfieldvalue(options,'vec_type','mpi'); +issmoptions.solver_type=getfieldvalue(options,'solver_type','mumps'); Index: /issm/trunk-jpl/src/m/solvers/issmsolver.py =================================================================== --- /issm/trunk-jpl/src/m/solvers/issmsolver.py (revision 16141) +++ /issm/trunk-jpl/src/m/solvers/issmsolver.py (revision 16142) @@ -9,5 +9,5 @@ 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. Index: /issm/trunk-jpl/test/NightlyRun/testad.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/testad.m (revision 16141) +++ /issm/trunk-jpl/test/NightlyRun/testad.m (revision 16142) @@ -1,11 +1,12 @@ %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;