Index: /issm/trunk/src/c/ConfigureObjectsx/ConfigureObjectsx.cpp =================================================================== --- /issm/trunk/src/c/ConfigureObjectsx/ConfigureObjectsx.cpp (revision 3372) +++ /issm/trunk/src/c/ConfigureObjectsx/ConfigureObjectsx.cpp (revision 3373) @@ -18,11 +18,14 @@ extern int my_rank; - //_printf_(" Configuring elements...\n"); + _printf_(" Configuring elements...\n"); elements->Configure(elements,loads,nodes,materials,parameters); - //_printf_(" Configuring loads...\n"); + _printf_(" Configuring loads...\n"); + //if (my_rank==2){ + // nodes->Echo(); + //} loads->Configure(elements,loads,nodes,materials,parameters); - //_printf_(" Configuring nodes...\n"); + _printf_(" Configuring nodes...\n"); nodes->Configure(elements,loads,nodes,materials,parameters); - //_printf_(" Configuring parameters...\n"); + _printf_(" Configuring parameters...\n"); parameters->Configure(elements,loads, nodes, materials,parameters); Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 3372) +++ /issm/trunk/src/c/Makefile.am (revision 3373) @@ -729,4 +729,5 @@ ./parallel/ProcessResults.cpp\ ./parallel/prognostic_core.cpp\ + ./parallel/prognostic2_core.cpp\ ./parallel/balancedthickness_core.cpp\ ./parallel/balancedvelocities_core.cpp\ @@ -804,5 +805,5 @@ bin_PROGRAMS = else -bin_PROGRAMS = diagnostic.exe thermal.exe prognostic.exe balancedthickness.exe balancedvelocities.exe transient.exe steadystate.exe slopecompute.exe +bin_PROGRAMS = diagnostic.exe thermal.exe prognostic.exe prognostic2.exe balancedthickness.exe balancedvelocities.exe transient.exe steadystate.exe slopecompute.exe endif @@ -821,4 +822,7 @@ prognostic_exe_CXXFLAGS= -fPIC -D_PARALLEL_ +prognostic2_exe_SOURCES = parallel/prognostic2.cpp +prognostic2_exe_CXXFLAGS= -fPIC -D_PARALLEL_ + balancedthickness_exe_SOURCES = parallel/balancedthickness.cpp balancedthickness_exe_CXXFLAGS= -fPIC -D_PARALLEL_ Index: /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp (revision 3372) +++ /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp (revision 3373) @@ -160,5 +160,5 @@ /*Used later on: */ - my_grids=(int*)xcalloc(iomodel->numberofnodes,sizeof(int)); + my_grids=(int*)xcalloc(3*iomodel->numberofelements,sizeof(int)); #endif @@ -223,7 +223,7 @@ into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids will hold which grids belong to this partition*/ - my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1; - my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1; - my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1; + my_grids[3*i+0]=1; + my_grids[3*i+1]=1; + my_grids[3*i+2]=1; #endif @@ -251,7 +251,7 @@ /*From the element partitioning, we can determine which grids are on the inside of this cpu's *element partition, and which are on its border with other nodes:*/ - gridborder=NewVec(iomodel->numberofnodes); - - for (i=0;inumberofnodes;i++){ + gridborder=NewVec(3*iomodel->numberofelements); + + for (i=0;i<3*iomodel->numberofelements;i++){ if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES); } Index: /issm/trunk/src/c/parallel/ProcessResults.cpp =================================================================== --- /issm/trunk/src/c/parallel/ProcessResults.cpp (revision 3372) +++ /issm/trunk/src/c/parallel/ProcessResults.cpp (revision 3373) @@ -123,4 +123,7 @@ if(analysis_type==PrognosticAnalysisEnum()){ fem_p=model->GetFormulation(PrognosticAnalysisEnum()); + } + if(analysis_type==Prognostic2AnalysisEnum()){ + fem_p=model->GetFormulation(Prognostic2AnalysisEnum()); } if(analysis_type==TransientAnalysisEnum()){ Index: /issm/trunk/src/c/parallel/parallel.h =================================================================== --- /issm/trunk/src/c/parallel/parallel.h (revision 3372) +++ /issm/trunk/src/c/parallel/parallel.h (revision 3373) @@ -17,4 +17,5 @@ void diagnostic_core(DataSet* results,Model* model, ParameterInputs* inputs); void prognostic_core(DataSet* results,Model* model, ParameterInputs* inputs); +void prognostic2_core(DataSet* results,Model* model, ParameterInputs* inputs); void balancedthickness_core(DataSet* results,Model* model, ParameterInputs* inputs); void balancedvelocities_core(DataSet* results,Model* model, ParameterInputs* inputs); Index: /issm/trunk/src/c/parallel/prognostic2.cpp =================================================================== --- /issm/trunk/src/c/parallel/prognostic2.cpp (revision 3373) +++ /issm/trunk/src/c/parallel/prognostic2.cpp (revision 3373) @@ -0,0 +1,165 @@ +/*!\file: prognostic2.cpp + * \brief: prognostic2 solution + */ + +#include "../issm.h" +#include "./parallel.h" + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + + +int main(int argc,char* *argv){ + + /*I/O: */ + FILE* fid=NULL; + char* inputfilename=NULL; + char* outputfilename=NULL; + char* lockname=NULL; + int numberofnodes; + double waitonlock=0; + + Model* model=NULL; + + double* vx_g=NULL; + double* vy_g=NULL; + double* h_g_initial=NULL; + double* m_g=NULL; + double* a_g=NULL; + double dt; + double yts; + int qmu_analysis; + + /*Results: */ + DataSet* results=NULL; + DataSet* processedresults=NULL; + Result* result=NULL; + + ParameterInputs* inputs=NULL; + Param* param=NULL; + + /*time*/ + double start, finish; + double start_core, finish_core; + double start_init, finish_init; + + MODULEBOOT(); + + #if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_)) + ISSMERROR(" parallel executable was compiled without support of parallel libraries!"); + #endif + + /*Initialize Petsc and get start time*/ + PetscInitialize(&argc,&argv,(char *)0,""); + MPI_Barrier(MPI_COMM_WORLD); start=MPI_Wtime(); + + /*Size and rank: */ + MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); + MPI_Comm_size(MPI_COMM_WORLD,&num_procs); + + _printf_("recover , input file name and output file name:\n"); + inputfilename=argv[2]; + outputfilename=argv[3]; + lockname=argv[4]; + + /*Initialize model structure: */ + MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime(); + model=new Model(); + + /*Open handle to data on disk: */ + fid=pfopen(inputfilename,"rb"); + + /*Initialize model structure: */ + model=new Model(); + + _printf_("read and create finite element model:\n"); + model->AddFormulation(fid,Prognostic2AnalysisEnum()); + + /*recover parameters: */ + model->FindParam(&waitonlock,"waitonlock"); + model->FindParam(&qmu_analysis,"qmu_analysis"); + + _printf_("initialize inputs:\n"); + + model->FindParam(&vx_g,NULL,NULL,"vx_g",Prognostic2AnalysisEnum()); + model->FindParam(&vy_g,NULL,NULL,"vy_g",Prognostic2AnalysisEnum()); + model->FindParam(&h_g_initial,NULL,NULL,"h_g",Prognostic2AnalysisEnum()); + model->FindParam(&m_g,NULL,NULL,"m_g",Prognostic2AnalysisEnum()); + model->FindParam(&a_g,NULL,NULL,"a_g",Prognostic2AnalysisEnum()); + model->FindParam(&dt,"dt"); + model->FindParam(&yts,"yts"); + model->FindParam(&numberofnodes,"numberofnodes"); + + inputs=new ParameterInputs; + inputs->Add("vx",vx_g,1,numberofnodes); + inputs->Add("vy",vy_g,1,numberofnodes); + inputs->Add("thickness",h_g_initial,1,numberofnodes); + inputs->Add("melting",m_g,1,numberofnodes); + inputs->Add("accumulation",a_g,1,numberofnodes); + inputs->Add("dt",dt*yts); + + _printf_("initialize results:\n"); + results=new DataSet(ResultsEnum()); + MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime(); + + /*are we running the solutoin sequence, or a qmu wrapper around it? : */ + if(!qmu_analysis){ + + /*run prognostic2 analysis: */ + _printf_("call computational core:\n"); + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + prognostic2_core(results,model,inputs); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + + } + else{ + + /*run qmu analysis: */ + _printf_("calling qmu analysis on prognostic2 core:\n"); + + #ifdef _HAVE_DAKOTA_ + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + Qmux(model,inputs,Prognostic2AnalysisEnum(),NoneAnalysisEnum()); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + #else + ISSMERROR(" Dakota not present, cannot do qmu!"); + #endif + } + + /*Add analysis_type to results: */ + result=new Result(results->Size()+1,0,1,"analysis_type","prognostic2"); + results->AddObject(result); + + _printf_("process results:\n"); + ProcessResults(&processedresults,results,model,Prognostic2AnalysisEnum()); + + _printf_("write results to disk:\n"); + OutputResults(processedresults,outputfilename); + + if (waitonlock>0){ + _printf_("write lock file:\n"); + WriteLockFile(lockname); + } + + /*Free ressources:*/ + delete processedresults; + delete results; + delete model; + delete inputs; + + /*Get finish time and close*/ + MPI_Barrier(MPI_COMM_WORLD); finish = MPI_Wtime( ); + _printf_("\n %-34s %f seconds \n","Model initialization elapsed time:",finish_init-start_init); + _printf_(" %-34s %f seconds \n","Core solution elapsed time:",finish_core-start_core); + _printf_("\n %s %i hrs %i min %i sec\n\n","Total elapsed time:",int((finish-start)/3600),int(int(finish-start)%3600/60),int(finish-start)%60); + _printf_("closing MPI and Petsc\n"); + PetscFinalize(); + + /*end module: */ + MODULEEND(); + + return 0; //unix success return; +} Index: /issm/trunk/src/c/parallel/prognostic2_core.cpp =================================================================== --- /issm/trunk/src/c/parallel/prognostic2_core.cpp (revision 3373) +++ /issm/trunk/src/c/parallel/prognostic2_core.cpp (revision 3373) @@ -0,0 +1,66 @@ +/*!\file: prognostic2_core.cpp + * \brief: core of the prognostic2 solution + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./parallel.h" +#include "../issm.h" + +void prognostic2_core(DataSet* results,Model* model,ParameterInputs* inputs){ + + extern int my_rank; + + /*output: */ + Result* result=NULL; + + /*intermediary: */ + Vec vx_g=NULL; + Vec vy_g=NULL; + + /*solutions: */ + Vec h_g=NULL; + + /*flags: */ + int verbose=0; + int numberofdofspernode; + int numberofnodes; + int dofs[1]={0}; + + /*fem prognostic model: */ + FemModel* fem_p=NULL; + + + /*recover fem model: */ + fem_p=model->GetFormulation(Prognostic2AnalysisEnum()); + + //first recover parameters common to all solutions + model->FindParam(&verbose,"verbose"); + model->FindParam(&numberofnodes,"numberofnodes"); + model->FindParam(&numberofdofspernode,"numberofdofspernode"); + + _printf_("depth averaging velocity...\n"); + vx_g=inputs->Get("vx",&dofs[0],1); + vy_g=inputs->Get("vy",&dofs[0],1); + FieldDepthAveragex(vx_g, fem_p->elements,fem_p->nodes, fem_p->loads, fem_p->materials,fem_p->parameters,"vx"); + FieldDepthAveragex(vy_g, fem_p->elements,fem_p->nodes, fem_p->loads, fem_p->materials,fem_p->parameters,"vy"); + inputs->Add("vx_average",vx_g,1,numberofnodes); + inputs->Add("vy_average",vy_g,1,numberofnodes); + + _printf_("call computational core:\n"); + diagnostic_core_linear(&h_g,fem_p,inputs,Prognostic2AnalysisEnum(),NoneAnalysisEnum()); + + _printf_("extrude computed thickness on all layers:\n"); + //FieldAverageOntoVertices(h_g, fem_p->elements,fem_p->nodes, fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0); + + /*Plug results into output dataset: */ + result=new Result(results->Size()+1,0,1,"h_g",h_g); + results->AddObject(result); + + /*Free ressources:*/ + VecFree(&vx_g); + VecFree(&vy_g); + VecFree(&h_g); +}