Index: /issm/trunk/src/c/solutions/ControlInitialization.cpp =================================================================== --- /issm/trunk/src/c/solutions/ControlInitialization.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/ControlInitialization.cpp (revision 3894) @@ -0,0 +1,175 @@ +/*!\file: ControlInitialization.cpp + * \brief: ... + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./parallel.h" +#include "../modules.h" + +void ControlInitialization(Model* model){ + + extern int my_rank; + + /*fem models: */ + FemModel* fem_dh=NULL; + FemModel* fem_dv=NULL; + FemModel* fem_dhu=NULL; + FemModel* fem_ds=NULL; + FemModel* fem_sl=NULL; + + /*solutions: */ + Vec ug=NULL; + Vec ug_horiz=NULL; + Vec ug_vert=NULL; + Vec ug_stokes=NULL; + Vec pg=NULL; + Vec slopex=NULL; + Vec slopey=NULL; + double* vx=NULL; + double* vy=NULL; + double* vz=NULL; + double* pressure=NULL; + + /*flags: */ + int verbose=0; + int dim=-1; + int ishutter=0; + int ismacayealpattyn=0; + int isstokes=0; + int numberofdofspernode_sl; + int numberofdofspernode_dh; + int numberofdofspernode_ds; + int numberofnodes; + + double stokesreconditioning; + + /*dof recovery: */ + int dof01[2]={0,1}; + int dof2[1]={2}; + int dof012[3]={0,1,2}; + int dof3[1]={3}; + double* dofset=NULL; + + /*first recover parameters common to all solutions:*/ + model->FindParam(&verbose,VerboseEnum); + model->FindParam(&dim,DimEnum); + model->FindParam(&ishutter,IsHutterEnum); + model->FindParam(&ismacayealpattyn,IsMacAyealPattynEnum); + model->FindParam(&isstokes,IsStokesEnum); + model->FindParam(&numberofnodes,NumberOfNodesEnum); + model->FindParam(&stokesreconditioning,StokesReconditioningEnum); + + /*recover fem models: */ + fem_dh=model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum); + fem_dv=model->GetFormulation(DiagnosticAnalysisEnum,VertAnalysisEnum); + fem_ds=model->GetFormulation(DiagnosticAnalysisEnum,StokesAnalysisEnum); + fem_dhu=model->GetFormulation(DiagnosticAnalysisEnum,HutterAnalysisEnum); + fem_sl=model->GetFormulation(SlopecomputeAnalysisEnum); + + //specific parameters for specific models + fem_dh->FindParam(&numberofdofspernode_dh,NumberOfDofsPerNodeEnum); + fem_sl->FindParam(&numberofdofspernode_sl,NumberOfDofsPerNodeEnum); + fem_ds->FindParam(&numberofdofspernode_ds,NumberOfDofsPerNodeEnum); + + /*if no Stokes, assign output and return*/ + if (!isstokes){ + model->SetActiveFormulation(fem_dh); + return; + } + + /*1: compute slopes once for all*/ + + //compute slopes + if(verbose)_printf_("%s\n","computing bed slope (x and y derivatives)..."); + diagnostic_core_linear(&slopex,fem_sl,SlopecomputeAnalysisEnum,BedXAnalysisEnum); + diagnostic_core_linear(&slopey,fem_sl,SlopecomputeAnalysisEnum,BedYAnalysisEnum); + FieldExtrudex( slopex, fem_sl->elements,fem_sl->nodes,fem_sl->vertices,fem_sl->loads,fem_sl->materials,fem_sl->parameters,"slopex",0); + FieldExtrudex( slopey, fem_sl->elements,fem_sl->nodes,fem_sl->vertices,fem_sl->loads,fem_sl->materials,fem_sl->parameters,"slopey",0); + + //Add in inputs + model->UpdateInputsFromVector(slopex,BedSlopexEnum,VertexEnum); + model->UpdateInputsFromVector(slopey,BedSlopeyEnum,VertexEnum); + VecFree(&slopex); VecFree(&slopey); + + /*2: run a complete diagnostic to update spcs*/ + + //horizontal velocity + if(verbose)_printf_("%s\n"," computing horizontal velocities..."); + diagnostic_core_nonlinear(&ug,NULL,NULL,NULL,fem_dh,DiagnosticAnalysisEnum,HorizAnalysisEnum); + if(verbose)_printf_("%s\n"," extruding horizontal velocities..."); + VecDuplicatePatch(&ug_horiz,ug); FieldExtrudex( ug_horiz,fem_dh->elements,fem_dh->nodes, fem_dh->vertices,fem_dh->loads,fem_dh-> materials,fem_dh->parameters,"velocity",1); + + //Add to inputs: + SplitSolutionVectorx(ug_horiz,numberofnodes,numberofdofspernode_dh,&vx,&vy); + model->UpdateInputsFromVector(vx,VxEnum,VertexEnum); + model->UpdateInputsFromVector(vy,VyEnum,VertexEnum); + + + //vertical velocity + if(verbose)_printf_("%s\n"," computing vertical velocities..."); + diagnostic_core_linear(&ug_vert,fem_dv,DiagnosticAnalysisEnum,VertAnalysisEnum); + + //Create 3d u_g + if(verbose)_printf_("%s\n"," combining horizontal and vertical velocities..."); + VecFree(&ug); ug=NewVec(numberofnodes*3); + xfree((void**)&dofset);dofset=dofsetgen(2,&dof01[0],3,numberofnodes*3); VecMerge(ug,ug_horiz,dofset,numberofnodes*2); + xfree((void**)&dofset);dofset=dofsetgen(1,&dof2[0],3,numberofnodes*3); VecMerge(ug,ug_vert,dofset,numberofnodes*1); + VecFree(&ug_vert); VecFree(&ug_horiz); + + //Create 4d u_g + if(verbose)_printf_("%s\n"," computing pressure according to Pattyn..."); + ComputePressurex( &pg,fem_dh->elements, fem_dh->nodes, fem_dh->vertices,fem_dh->loads, fem_dh->materials, fem_dh->parameters,DiagnosticAnalysisEnum,HorizAnalysisEnum); + VecScale(pg,1.0/stokesreconditioning); + ug_stokes=NewVec(fem_ds->nodesets->GetGSize()); + xfree((void**)&dofset);dofset=dofsetgen(3,dof012,4,numberofnodes*4); VecMerge(ug_stokes,ug,dofset,numberofnodes*3); + xfree((void**)&dofset);dofset=dofsetgen(1,dof3,4,numberofnodes*4); VecMerge(ug_stokes,pg,dofset,numberofnodes); + + //Add in inputs + SplitSolutionVectorx(ug_stokes,numberofnodes,numberofdofspernode_ds,&vx,&vy,&vz,&pressure); + model->UpdateInputsFromVector(vx,VxEnum,VertexEnum); + model->UpdateInputsFromVector(vy,VyEnum,VertexEnum); + model->UpdateInputsFromVector(vz,VzEnum,VertexEnum); + model->UpdateInputsFromVector(pressure,PressureEnum,VertexEnum); + VecFree(&ug_stokes); + + //update spcs + if(verbose)_printf_("%s\n"," update boundary conditions for stokes using velocities previously computed..."); + xfree((void**)&dofset);dofset=dofsetgen(3,dof012,4,numberofnodes*4); VecMerge(fem_ds->yg->vector,ug,dofset,3*numberofnodes); + + VecFree(&fem_ds->ys); VecFree(&fem_ds->ys0); + Reducevectorgtosx(&fem_ds->ys,&fem_ds->ys0, fem_ds->yg->vector,fem_ds->nodesets); + + //Compute Stokes velocities to speed up later runs + if(verbose)_printf_("%s\n"," computing stokes velocities and pressure ..."); + VecFree(&ug); + diagnostic_core_nonlinear(&ug,NULL,NULL,NULL,fem_ds,DiagnosticAnalysisEnum,StokesAnalysisEnum); + + //Add in inputs + xfree((void**)&vx); xfree((void**)&vy); xfree((void**)&vz); xfree((void**)&pressure); + SplitSolutionVectorx(ug,numberofnodes,numberofdofspernode_ds,&vx,&vy,&vz,&pressure); + model->UpdateInputsFromVector(vx,VxEnum,VertexEnum); + model->UpdateInputsFromVector(vy,VyEnum,VertexEnum); + model->UpdateInputsFromVector(vz,VzEnum,VertexEnum); + model->UpdateInputsFromVector(pressure,PressureEnum,VertexEnum); + + /*Assign output*/ + model->SetActiveFormulation(fem_ds); + + /*Free ressources:*/ + xfree((void**)&dofset); + xfree((void**)&vx); + xfree((void**)&vy); + xfree((void**)&vz); + xfree((void**)&pressure); + VecFree(&ug); + VecFree(&ug_horiz); + VecFree(&ug_vert); + VecFree(&ug_stokes); + VecFree(&pg); + VecFree(&slopex); + VecFree(&slopey); + +} Index: /issm/trunk/src/c/solutions/ControlRestart.cpp =================================================================== --- /issm/trunk/src/c/solutions/ControlRestart.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/ControlRestart.cpp (revision 3894) @@ -0,0 +1,55 @@ + +/*!\file: ControlRestart.cpp + * \brief: core of the control solution + */ + +#include "./parallel.h" +#include "../modules.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +void ControlRestart(Model* model,double* param_g){ + + extern int my_rank; + + /*output: */ + DataSet* temporary_results=NULL; + DataSet* results=NULL; + Result* result=NULL; + char* outputfilename=NULL; + + /*Intermediary: */ + int i; + int numberofnodes; + double* param_g_copy; + + /*Recover parameters used throughout the solution:*/ + model->FindParam(&numberofnodes,NumberOfNodesEnum); + model->FindParam(&outputfilename,OutputFileNameEnum); + + /*Plug COPYS of the results into output dataset: + * only the pointer is given to temporary_results and at the + * end of ProcessResults the pointer is deleted. That would + * destroy param_g*/ + + param_g_copy=(double*)xcalloc(numberofnodes,sizeof(double)); + for(i=0;iSize()+1,0,1,"param_g",param_g_copy,numberofnodes); + temporary_results->AddObject(result); + + result=new Result(temporary_results->Size()+1,0,1,"analysis_type",EnumAsString(DiagnosticAnalysisEnum)); + temporary_results->AddObject(result); + + //process results + ProcessResults(&results,temporary_results,model,ControlAnalysisEnum); + + //Write results on disk + OutputResults(results,outputfilename); + + /*Free ressources:*/ + delete temporary_results; + delete results; + xfree((void**)&outputfilename); + xfree((void**)¶m_g_copy); +} Index: /issm/trunk/src/c/solutions/OutputResults.cpp =================================================================== --- /issm/trunk/src/c/solutions/OutputResults.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/OutputResults.cpp (revision 3894) @@ -0,0 +1,37 @@ +/*!\file: OutputResults.cpp + * \brief: go through results dataset, and for each result, write it to disk. + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "stdio.h" +#include "../DataSet/DataSet.h" +#include "../io/io.h" +#include "../objects/objects.h" + +void OutputResults(DataSet* results,char* filename){ + + int i; + Result* result=NULL; + FILE* fid=NULL; + extern int my_rank; + + /* Open output file to write raw binary data: */ + if(my_rank==0){ + fid=pfopen(filename,"wb"); + + for(i=0;iSize();i++){ + result=(Result*)results->GetObjectByOffset(i); + + /*write result to disk: */ + result->WriteData(fid); + } + + /*Close file: */ + pfclose(fid,filename); + } +} Index: /issm/trunk/src/c/solutions/ProcessResults.cpp =================================================================== --- /issm/trunk/src/c/solutions/ProcessResults.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/ProcessResults.cpp (revision 3894) @@ -0,0 +1,546 @@ +/*!\file: ProcessResults.cpp + * \brief: go through results dataset, and for each result, process it for easier retrieval + * by the Matlab side. This usually means splitting the velocities from the g-size nodeset + * to the grid set (ug->vx,vy,vz), same for pressure (p_g->pressure), etc ... It also implies + * departitioning of the results. + * This phase is necessary prior to outputting the results on disk. + */ + +#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 "../objects/objects.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../shared/shared.h" + +void ProcessResults(DataSet** pnewresults, DataSet* results,Model* model,int analysis_type){ + + int i,n; + Result* result=NULL; + Result* newresult=NULL; + + /*output: */ + DataSet* newresults=NULL; + + /*fem diagnostic models: */ + FemModel* fem_dh=NULL; + FemModel* fem_dv=NULL; + FemModel* fem_dhu=NULL; + FemModel* fem_ds=NULL; + FemModel* fem_sl=NULL; + + /*fem thermal models: */ + FemModel* fem_t=NULL; + + /*fem prognostic models: */ + FemModel* fem_p=NULL; + + /*fem control models: */ + FemModel* fem_c=NULL; + + /*some parameters*/ + bool ishutter; + bool ismacayealpattyn; + bool isstokes; + int dim; + + /*intermediary: */ + Vec u_g=NULL; + double* u_g_serial=NULL; + double* vx=NULL; + double* vy=NULL; + double* vz=NULL; + double* vel=NULL; + Vec p_g=NULL; + double* p_g_serial=NULL; + double* pressure=NULL; + double* partition=NULL; + double yts; + + double* param_g=NULL; + double* parameter=NULL; + Vec riftproperties=NULL; + double* riftproperties_serial=NULL; + int numrifts=0; + + Vec s_g=NULL; + double* s_g_serial=NULL; + double* surface=NULL; + + Vec sx_g=NULL; + double* sx_g_serial=NULL; + double* slopex=NULL; + + Vec sy_g=NULL; + double* sy_g_serial=NULL; + double* slopey=NULL; + + Vec b_g=NULL; + double* b_g_serial=NULL; + double* bed=NULL; + + Vec h_g=NULL; + double* h_g_serial=NULL; + double* thickness=NULL; + + Vec t_g=NULL; + double* t_g_serial=NULL; + double* temperature=NULL; + + Vec m_g=NULL; + double* m_g_serial=NULL; + double* melting=NULL; + + Vec grad_g=NULL; + double* grad_g_serial=NULL; + double* gradient=NULL; + + Vec sigma_zz=NULL; + double* sigma_zz_serial=NULL; + + Vec v_g=NULL; + double* v_g_serial=NULL; + + int numberofnodes,numberofvertices,numberofelements; + + /*Initialize new results: */ + newresults=new DataSet(ResultsEnum); + + /*some flags needed: */ + model->FindParam(&dim,DimEnum); + model->FindParam(&ishutter,IsHutterEnum); + model->FindParam(&isstokes,IsStokesEnum); + model->FindParam(&ismacayealpattyn,IsMacAyealPattynEnum); + + /*Recover femmodels first: */ + fem_dh=model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum); + fem_c=fem_dh; + fem_dv=model->GetFormulation(DiagnosticAnalysisEnum,VertAnalysisEnum); + fem_ds=model->GetFormulation(DiagnosticAnalysisEnum,StokesAnalysisEnum); + fem_dhu=model->GetFormulation(DiagnosticAnalysisEnum,HutterAnalysisEnum); + fem_sl=model->GetFormulation(SlopecomputeAnalysisEnum); + if(analysis_type==PrognosticAnalysisEnum){ + fem_p=model->GetFormulation(PrognosticAnalysisEnum); + } + if(analysis_type==Prognostic2AnalysisEnum){ + fem_p=model->GetFormulation(Prognostic2AnalysisEnum); + } + if(analysis_type==TransientAnalysisEnum){ + fem_p=model->GetFormulation(PrognosticAnalysisEnum); + } + if(analysis_type==BalancedthicknessAnalysisEnum){ + fem_p=model->GetFormulation(BalancedthicknessAnalysisEnum); + } + if(analysis_type==Balancedthickness2AnalysisEnum){ + fem_p=model->GetFormulation(Balancedthickness2AnalysisEnum); + } + if(analysis_type==BalancedvelocitiesAnalysisEnum){ + fem_p=model->GetFormulation(BalancedvelocitiesAnalysisEnum); + } + fem_t=model->GetFormulation(ThermalAnalysisEnum); + + for(n=0;nSize();n++){ + result=(Result*)results->GetObjectByOffset(n); + + if(strcmp(result->GetFieldName(),"u_g")==0){ + + /*Ok, are we dealing with velocities coming from MacAyeal, Pattyin, Hutter, on 2,3 dofs or + *Stokes on 4 dofs: */ + result->GetField(&u_g); + VecToMPISerial(&u_g_serial,u_g); + + //2d results -> 2 dofs per node + if (dim==2){ + /*ok, 2 dofs, on number of nodes: */ + if(ismacayealpattyn){ + fem_dh->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + VecToMPISerial(&partition,fem_dh->partition->vector); + fem_dh->parameters->FindParam(&yts,YtsEnum); + } + else{ + fem_dhu->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + VecToMPISerial(&partition,fem_dhu->partition->vector); + fem_dhu->parameters->FindParam(&yts,YtsEnum); + } + vx=(double*)xmalloc(numberofnodes*sizeof(double)); + vy=(double*)xmalloc(numberofnodes*sizeof(double)); + vz=(double*)xmalloc(numberofnodes*sizeof(double)); + vel=(double*)xmalloc(numberofnodes*sizeof(double)); + + for(i=0;i 3 or 4 (stokes) dofs per node + else{ + if(!isstokes){ + /*ok, 3 dofs, on number of nodes: */ + if(ismacayealpattyn){ + fem_dh->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + VecToMPISerial(&partition,fem_dh->partition->vector); + fem_dh->parameters->FindParam(&yts,YtsEnum); + } + else{ + fem_dhu->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + VecToMPISerial(&partition,fem_dhu->partition->vector); + fem_dhu->parameters->FindParam(&yts,YtsEnum); + } + vx=(double*)xmalloc(numberofnodes*sizeof(double)); + vy=(double*)xmalloc(numberofnodes*sizeof(double)); + vz=(double*)xmalloc(numberofnodes*sizeof(double)); + vel=(double*)xmalloc(numberofnodes*sizeof(double)); + + for(i=0;iparameters->FindParam(&numberofnodes,NumberOfNodesEnum); + VecToMPISerial(&partition,fem_ds->partition->vector); + fem_ds->parameters->FindParam(&yts,YtsEnum); + vx=(double*)xmalloc(numberofnodes*sizeof(double)); + vy=(double*)xmalloc(numberofnodes*sizeof(double)); + vz=(double*)xmalloc(numberofnodes*sizeof(double)); + vel=(double*)xmalloc(numberofnodes*sizeof(double)); + for(i=0;iSize()+1,result->GetTime(),result->GetStep(),"vx",vx,numberofnodes); + newresults->AddObject(newresult); + + newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vy",vy,numberofnodes); + newresults->AddObject(newresult); + + newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vz",vz,numberofnodes); + 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); + xfree((void**)&partition); + xfree((void**)&vx); + xfree((void**)&vy); + xfree((void**)&vz); + xfree((void**)&vel); + VecFree(&u_g); + } + else if(strcmp(result->GetFieldName(),"p_g")==0){ + /*easy, p_g is of size numberofnodes, on 1 dof, just repartition: */ + result->GetField(&p_g); + VecToMPISerial(&p_g_serial,p_g); + + if(!isstokes){ + if(ismacayealpattyn){ + fem_dh->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + VecToMPISerial(&partition,fem_dh->partition->vector); + } + else{ + fem_dhu->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + VecToMPISerial(&partition,fem_dhu->partition->vector); + } + } + else{ + fem_ds->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + VecToMPISerial(&partition,fem_ds->partition->vector); + } + + pressure=(double*)xmalloc(numberofnodes*sizeof(double)); + + for(i=0;iSize()+1,result->GetTime(),result->GetStep(),"pressure",pressure,numberofnodes); + newresults->AddObject(newresult); + + /*do some cleanup: */ + xfree((void**)&p_g_serial); + xfree((void**)&partition); + xfree((void**)&pressure); + VecFree(&p_g); + } + else if(strcmp(result->GetFieldName(),"t_g")==0){ + /*easy, t_g is of size numberofnodes, on 1 dof, just repartition: */ + result->GetField(&t_g); + VecToMPISerial(&t_g_serial,t_g); + fem_t->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + VecToMPISerial(&partition,fem_t->partition->vector); + + temperature=(double*)xmalloc(numberofnodes*sizeof(double)); + + for(i=0;iSize()+1,result->GetTime(),result->GetStep(),"temperature",temperature,numberofnodes); + newresults->AddObject(newresult); + + /*do some cleanup: */ + xfree((void**)&t_g_serial); + xfree((void**)&partition); + xfree((void**)&temperature); + VecFree(&t_g); + } + else if(strcmp(result->GetFieldName(),"grad_g")==0){ + + /*easy, grad_g is of size numberofnodes, on 1 dof, just repartition: */ + result->GetField(&grad_g); + VecToMPISerial(&grad_g_serial,grad_g); + fem_c->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + VecToMPISerial(&partition,fem_c->partition->vector); + + gradient=(double*)xmalloc(numberofnodes*sizeof(double)); + + for(i=0;iSize()+1,result->GetTime(),result->GetStep(),"gradient",gradient,numberofnodes); + newresults->AddObject(newresult); + + /*do some cleanup: */ + xfree((void**)&grad_g_serial); + xfree((void**)&partition); + xfree((void**)&gradient); + VecFree(&grad_g); + } + else if(strcmp(result->GetFieldName(),"m_g")==0){ + /*easy, m_g is of size numberofnodes, on 1 dof, just repartition: */ + result->GetField(&m_g); + VecToMPISerial(&m_g_serial,m_g); + fem_t->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + fem_t->parameters->FindParam(&yts,YtsEnum); + VecToMPISerial(&partition,fem_t->partition->vector); + + melting=(double*)xmalloc(numberofnodes*sizeof(double)); + + for(i=0;iSize()+1,result->GetTime(),result->GetStep(),"melting",melting,numberofnodes); + newresults->AddObject(newresult); + + /*do some cleanup: */ + xfree((void**)&m_g_serial); + xfree((void**)&partition); + xfree((void**)&melting); + VecFree(&m_g); + } + else if(strcmp(result->GetFieldName(),"h_g")==0){ + /*easy, h_g is of size numberofvertices, on 1 dof, just repartition: */ + result->GetField(&h_g); + VecToMPISerial(&h_g_serial,h_g); + fem_p->parameters->FindParam(&numberofvertices,NumberOfVerticesEnum); + VecToMPISerial(&partition,fem_p->partition->vector); + + thickness=(double*)xmalloc(numberofvertices*sizeof(double)); + + for(i=0;iSize()+1,result->GetTime(),result->GetStep(),"thickness",thickness,numberofvertices); + newresults->AddObject(newresult); + + /*do some cleanup: */ + xfree((void**)&h_g_serial); + xfree((void**)&thickness); + xfree((void**)&partition); + VecFree(&h_g); + } + else if(strcmp(result->GetFieldName(),"v_g")==0){ + /*easy, v_g is of size numberofnodes, on 1 dof, just repartition: */ + result->GetField(&v_g); + VecToMPISerial(&v_g_serial,v_g); + fem_p->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + VecToMPISerial(&partition,fem_p->partition->vector); + + vel=(double*)xmalloc(numberofnodes*sizeof(double)); + + for(i=0;iSize()+1,result->GetTime(),result->GetStep(),"vel",vel,numberofnodes); + newresults->AddObject(newresult); + + /*do some cleanup: */ + xfree((void**)&v_g_serial); + xfree((void**)&vel); + xfree((void**)&partition); + VecFree(&v_g); + } + else if(strcmp(result->GetFieldName(),"s_g")==0){ + /*easy, s_g is of size numberofnodes, on 1 dof, just repartition: */ + result->GetField(&s_g); + VecToMPISerial(&s_g_serial,s_g); + fem_p->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + VecToMPISerial(&partition,fem_p->partition->vector); + + surface=(double*)xmalloc(numberofnodes*sizeof(double)); + + for(i=0;iSize()+1,result->GetTime(),result->GetStep(),"surface",surface,numberofnodes); + newresults->AddObject(newresult); + + /*do some cleanup: */ + xfree((void**)&s_g_serial); + xfree((void**)&partition); + xfree((void**)&surface); + VecFree(&s_g); + } + else if(strcmp(result->GetFieldName(),"sx_g")==0){ + /*easy, s_g is of size numberofnodes, on 1 dof, just repartition: */ + result->GetField(&sx_g); + VecToMPISerial(&sx_g_serial,sx_g); + fem_sl->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + VecToMPISerial(&partition,fem_sl->partition->vector); + + slopex=(double*)xmalloc(numberofnodes*sizeof(double)); + + for(i=0;iSize()+1,result->GetTime(),result->GetStep(),"slopex",slopex,numberofnodes); + newresults->AddObject(newresult); + + /*do some cleanup: */ + xfree((void**)&sx_g_serial); + xfree((void**)&partition); + xfree((void**)&slopex); + VecFree(&sx_g); + } + else if(strcmp(result->GetFieldName(),"sy_g")==0){ + /*easy, s_g is of size numberofnodes, on 1 dof, just repartition: */ + result->GetField(&sy_g); + VecToMPISerial(&sy_g_serial,sy_g); + fem_sl->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + VecToMPISerial(&partition,fem_sl->partition->vector); + + slopey=(double*)xmalloc(numberofnodes*sizeof(double)); + + for(i=0;iSize()+1,result->GetTime(),result->GetStep(),"slopey",slopey,numberofnodes); + newresults->AddObject(newresult); + + /*do some cleanup: */ + xfree((void**)&sy_g_serial); + xfree((void**)&partition); + xfree((void**)&slopey); + VecFree(&sy_g); + } + else if(strcmp(result->GetFieldName(),"b_g")==0){ + /*easy, b_g is of size numberofnodes, on 1 dof, just repartition: */ + result->GetField(&b_g); + VecToMPISerial(&b_g_serial,b_g); + fem_p->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + VecToMPISerial(&partition,fem_p->partition->vector); + + bed=(double*)xmalloc(numberofnodes*sizeof(double)); + + for(i=0;iSize()+1,result->GetTime(),result->GetStep(),"bed",bed,numberofnodes); + newresults->AddObject(newresult); + + /*do some cleanup: */ + xfree((void**)&b_g_serial); + xfree((void**)&partition); + xfree((void**)&bed); + VecFree(&b_g); + } + else if(strcmp(result->GetFieldName(),"param_g")==0){ + /*easy, param_g is of size numberofnodes, on 1 dof, just repartition: */ + result->GetField(¶m_g); + fem_dh->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + VecToMPISerial(&partition,fem_dh->partition->vector); + + parameter=(double*)xmalloc(numberofnodes*sizeof(double)); + + for(i=0;iSize()+1,result->GetTime(),result->GetStep(),"parameter",parameter,numberofnodes); + newresults->AddObject(newresult); + + /*do some cleanup: */ + xfree((void**)&partition); + xfree((void**)¶m_g); + xfree((void**)¶meter); + } + else if(strcmp(result->GetFieldName(),"riftproperties")==0){ + result->GetField(&riftproperties); + fem_dh->parameters->FindParam(&numrifts,NumRiftsEnum); + VecToMPISerial(&riftproperties_serial,riftproperties); + + /*Ok, add parameter to newresults: */ + newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"riftproperties",riftproperties_serial,numrifts); + newresults->AddObject(newresult); + xfree((void**)&riftproperties); + + } + else if(strcmp(result->GetFieldName(),"sigma_zz")==0){ + /*easy, param_g is of size numberofelements, on 1 dof, just repartition: */ + fem_ds->parameters->FindParam(&numberofelements,NumberOfElementsEnum); + result->GetField(&sigma_zz); + VecToMPISerial(&sigma_zz_serial,sigma_zz); + + /*Ok, add parameter to newresults: */ + newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"sigma_zz",sigma_zz_serial,numberofelements); + newresults->AddObject(newresult); + + /*do some cleanup: */ + xfree((void**)&sigma_zz_serial); + + } + else{ + /*Just copy the result into the new results dataset: */ + newresults->AddObject(result->copy()); + } + } + + /*Assign output pointers:*/ + *pnewresults=newresults; +} Index: /issm/trunk/src/c/solutions/WriteLockFile.cpp =================================================================== --- /issm/trunk/src/c/solutions/WriteLockFile.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/WriteLockFile.cpp (revision 3894) @@ -0,0 +1,25 @@ +/*!\file: WriteLockFile.cpp + * \brief + */ +#include "stdio.h" +#include "../shared/shared.h" +#include "../include/include.h" + +void WriteLockFile(char* filename){ + + int i; + extern int my_rank; + + /* output: */ + FILE* fid=NULL; + + /* Open lock file and write 1 into it: */ + if(my_rank==0){ + fid=fopen(filename,"w"); + if(fid==NULL) ISSMERROR("%s%s","error message: could not open lock file ",filename); + + /*Close file: */ + if(fclose(fid)!=0) ISSMERROR("%s%s","could not close lock file ",filename); + } + +} Index: /issm/trunk/src/c/solutions/balancedthickness.cpp =================================================================== --- /issm/trunk/src/c/solutions/balancedthickness.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/balancedthickness.cpp (revision 3894) @@ -0,0 +1,135 @@ +/*!\file: balancedthickness.cpp + * \brief: balancedthickness solution + */ +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../DataSet/DataSet.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../include/include.h" +#include "../modules.h" +#include "./parallel.h" + +int main(int argc,char* *argv){ + + /*I/O: */ + FILE* fid=NULL; + char* inputfilename=NULL; + char* outputfilename=NULL; + char* lockname=NULL; + int numberofnodes; + bool waitonlock=false; + + Model* model=NULL; + + double dt; + double yts; + bool qmu_analysis; + + /*Results: */ + DataSet* results=NULL; + DataSet* processedresults=NULL; + Result* result=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(); + + /*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,BalancedthicknessAnalysisEnum); + + /*recover parameters: */ + model->FindParam(&waitonlock,WaitOnLockEnum); + model->FindParam(&qmu_analysis,QmuAnalysisEnum); + + 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 balancedthickness analysis: */ + _printf_("call computational core:\n"); + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + results=balancedthickness_core(model); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + + } + else{ + + /*run qmu analysis: */ + _printf_("calling qmu analysis on balancedthickness core:\n"); + + #ifdef _HAVE_DAKOTA_ + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + Qmux(model,BalancedthicknessAnalysisEnum,NoneAnalysisEnum); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + #else + ISSMERROR(" Dakota not present, cannot do qmu!"); + #endif + } + + _printf_("process results:\n"); + ProcessResults(&processedresults,results,model,BalancedthicknessAnalysisEnum); + + _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; + + /*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/solutions/balancedthickness2.cpp =================================================================== --- /issm/trunk/src/c/solutions/balancedthickness2.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/balancedthickness2.cpp (revision 3894) @@ -0,0 +1,135 @@ +/*!\file: balancedthickness2.cpp + * \brief: balancedthickness2 solution + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../DataSet/DataSet.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../include/include.h" +#include "../modules.h" +#include "./parallel.h" + +int main(int argc,char* *argv){ + + /*I/O: */ + FILE* fid=NULL; + char* inputfilename=NULL; + char* outputfilename=NULL; + char* lockname=NULL; + bool waitonlock=false; + + Model* model=NULL; + + double dt; + double yts; + bool qmu_analysis; + + /*Results: */ + DataSet* results=NULL; + DataSet* processedresults=NULL; + Result* result=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(); + + /*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,Balancedthickness2AnalysisEnum); + + /*recover parameters: */ + model->FindParam(&waitonlock,WaitOnLockEnum); + model->FindParam(&qmu_analysis,QmuAnalysisEnum); + + 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 balancedthickness2 analysis: */ + _printf_("call computational core:\n"); + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + results=balancedthickness2_core(model); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + + } + else{ + + /*run qmu analysis: */ + _printf_("calling qmu analysis on balancedthickness2 core:\n"); + + #ifdef _HAVE_DAKOTA_ + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + Qmux(model,Balancedthickness2AnalysisEnum,NoneAnalysisEnum); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + #else + ISSMERROR(" Dakota not present, cannot do qmu!"); + #endif + } + + _printf_("process results:\n"); + ProcessResults(&processedresults,results,model,Balancedthickness2AnalysisEnum); + + _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; + + /*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/solutions/balancedthickness2_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/balancedthickness2_core.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/balancedthickness2_core.cpp (revision 3894) @@ -0,0 +1,85 @@ +/*!\file: balancedthickness2_core.cpp + * \brief: core of the balancedthickness2 solution + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../include/include.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./parallel.h" +#include "../modules.h" + +DataSet* balancedthickness2_core(Model* model){ + + extern int my_rank; + + /*output: */ + Result* result=NULL; + DataSet* results=NULL; + + /*intermediary: */ + Vec vx_g=NULL; + Vec vy_g=NULL; + + /*solutions: */ + Vec h_g=NULL; + + /*flags: */ + int verbose=0; + int numberofdofspernode; + int numberofnodes; + int numberofvertices; + int dofs[1]={1}; + + /*fem balancedthickness2 model: */ + FemModel* fem_p=NULL; + + //initialize results: + results=new DataSet(ResultsEnum); + + fem_p=model->GetFormulation(Balancedthickness2AnalysisEnum); + + //first recover parameters common to all solutions + model->FindParam(&verbose,VerboseEnum); + model->FindParam(&numberofnodes,NumberOfNodesEnum); + model->FindParam(&numberofvertices,NumberOfVerticesEnum); + model->FindParam(&numberofdofspernode,NumberOfDofsPerNodeEnum); + + _printf_("depth averaging velocity...\n"); + //vx_g=inputs->Get("vx",&dofs[0],1); + //vy_g=inputs->Get("vy",&dofs[0],1); + /* NOT WORKING YET.... + FieldDepthAveragex(vx_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"vx"); + FieldDepthAveragex(vy_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"vy"); + */ + //inputs->Add("vx_average",vx_g,1,numberofvertices); + //inputs->Add("vy_average",vy_g,1,numberofvertices); + ISSMERROR("not supported yet!"); + + _printf_("call computational core:\n"); + diagnostic_core_linear(&h_g,fem_p,Balancedthickness2AnalysisEnum,NoneAnalysisEnum); + + _printf_("Averaging over vertices:\n"); + FieldAverageOntoVerticesx(&h_g,fem_p->elements,fem_p->nodes,fem_p->vertices,fem_p->loads,fem_p->materials,fem_p->parameters); + +// _printf_("extrude computed thickness on all layers:\n"); +// FieldExtrudex( h_g, fem_p->elements,fem_p->nodes, fem_p->vertices,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); + + /*Add analysis_type to results: */ + result=new Result(results->Size()+1,0,1,"analysis_type","balancedthickness2"); + results->AddObject(result); + + /*Free ressources:*/ + VecFree(&vx_g); + VecFree(&vy_g); + VecFree(&h_g); + + /*return: */ + return results; + +} Index: /issm/trunk/src/c/solutions/balancedthickness_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/balancedthickness_core.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/balancedthickness_core.cpp (revision 3894) @@ -0,0 +1,74 @@ +/*!\file: balancedthickness_core.cpp + * \brief: core of the balancedthickness solution + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../include/include.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./parallel.h" +#include "../modules.h" + +DataSet* balancedthickness_core(Model* model){ + + extern int my_rank; + + /*output: */ + Result* result=NULL; + DataSet* results=NULL; + + /*intermediary: */ + Vec u_g=NULL; + + /*solutions: */ + Vec h_g=NULL; + + /*flags: */ + int verbose=0; + int numberofdofspernode; + int numberofnodes; + int dofs[2]={1,1}; + + /*fem balancedthickness model: */ + FemModel* fem_p=NULL; + + //initialize results: + results=new DataSet(ResultsEnum); + + /*recover fem model: */ + fem_p=model->GetFormulation(BalancedthicknessAnalysisEnum); + + //first recover parameters common to all solutions + model->FindParam(&verbose,VerboseEnum); + model->FindParam(&numberofnodes,NumberOfNodesEnum); + model->FindParam(&numberofdofspernode,NumberOfDofsPerNodeEnum); + + _printf_("depth averaging velocity...\n"); + //u_g=inputs->Get("velocity",&dofs[0],2); //take (vx,vy) from inputs velocity + //FieldDepthAveragex( u_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"velocity"); + //inputs->Add("velocity_average",u_g,2,numberofnodes); + ISSMERROR("not supported yet!"); + + _printf_("call computational core:\n"); + diagnostic_core_linear(&h_g,fem_p,BalancedthicknessAnalysisEnum,NoneAnalysisEnum); + + _printf_("extrude computed thickness on all layers:\n"); + FieldExtrudex( h_g, fem_p->elements,fem_p->nodes, fem_p->vertices,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); + + /*Add analysis_type to results: */ + result=new Result(results->Size()+1,0,1,"analysis_type","balancedthickness"); + results->AddObject(result); + + /*Free ressources:*/ + VecFree(&u_g); + VecFree(&h_g); + + /*return: */ + return results; + +} Index: /issm/trunk/src/c/solutions/balancedvelocities.cpp =================================================================== --- /issm/trunk/src/c/solutions/balancedvelocities.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/balancedvelocities.cpp (revision 3894) @@ -0,0 +1,133 @@ +/*!\file: balancedvelocities.cpp + * \brief: balancedvelocities solution + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../DataSet/DataSet.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../include/include.h" +#include "../modules.h" +#include "./parallel.h" + +int main(int argc,char* *argv){ + + /*I/O: */ + FILE* fid=NULL; + char* inputfilename=NULL; + char* outputfilename=NULL; + char* lockname=NULL; + int numberofnodes; + bool waitonlock=false; + + Model* model=NULL; + + bool qmu_analysis; + + /*Results: */ + DataSet* results=NULL; + DataSet* processedresults=NULL; + Result* result=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(); + + /*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,BalancedvelocitiesAnalysisEnum); + + /*recover parameters: */ + model->FindParam(&waitonlock,WaitOnLockEnum); + model->FindParam(&qmu_analysis,QmuAnalysisEnum); + + 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 balancedvelocities analysis: */ + _printf_("call computational core:\n"); + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + results=balancedvelocities_core(model); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + + } + else{ + + /*run qmu analysis: */ + _printf_("calling qmu analysis on balancedvelocities core:\n"); + + #ifdef _HAVE_DAKOTA_ + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + Qmux(model,BalancedvelocitiesAnalysisEnum,NoneAnalysisEnum); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + #else + ISSMERROR(" Dakota not present, cannot do qmu!"); + #endif + } + + _printf_("process results:\n"); + ProcessResults(&processedresults,results,model,BalancedvelocitiesAnalysisEnum); + + _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; + + /*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/solutions/balancedvelocities_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/balancedvelocities_core.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/balancedvelocities_core.cpp (revision 3894) @@ -0,0 +1,72 @@ +/*!\file: balancedvelocities_core.cpp + * \brief: core of the balancedvelocities solution + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./parallel.h" +#include "../modules.h" + +DataSet* balancedvelocities_core(Model* model){ + + extern int my_rank; + + /*output: */ + Result* result=NULL; + DataSet* results=NULL; + + /*intermediary: */ + Vec u_g=NULL; + + /*solutions: */ + Vec v_g=NULL; + + /*flags: */ + int verbose=0; + int numberofdofspernode; + int numberofnodes; + int dofs[2]={1,1}; + + /*fem balancedvelocities model: */ + FemModel* fem_p=NULL; + + //initialize results: + results=new DataSet(ResultsEnum); + + /*recover fem model: */ + fem_p=model->GetFormulation(BalancedvelocitiesAnalysisEnum); + + //first recover parameters common to all solutions + model->FindParam(&verbose,VerboseEnum); + model->FindParam(&numberofnodes,NumberOfNodesEnum); + model->FindParam(&numberofdofspernode,NumberOfDofsPerNodeEnum); + + _printf_("depth averaging velocity...\n"); + //u_g=inputs->Get("velocity",&dofs[0],2); //take (vx,vy) from inputs velocity + //FieldDepthAveragex( u_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"velocity"); + //inputs->Add("velocity_average",u_g,2,numberofnodes); + ISSMERROR(" not supported yet!"); + + _printf_("call computational core:\n"); + diagnostic_core_linear(&v_g,fem_p,BalancedvelocitiesAnalysisEnum,NoneAnalysisEnum); + + _printf_("extrude computed thickness on all layers:\n"); + FieldExtrudex( v_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0); + + /*Plug results into output dataset: */ + result=new Result(results->Size()+1,0,1,"v_g",v_g); + results->AddObject(result); + + /*Add analysis_type to results: */ + result=new Result(results->Size()+1,0,1,"analysis_type","balancedvelocities"); + results->AddObject(result); + + /*Free ressources:*/ + VecFree(&u_g); + VecFree(&v_g); + + /*return: */ + return results; +} Index: /issm/trunk/src/c/solutions/control_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/control_core.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/control_core.cpp (revision 3894) @@ -0,0 +1,257 @@ +/*!\file: control_core.cpp + * \brief: core of the control solution + */ + +#include "./parallel.h" +#include "../modules.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +DataSet* control_core(Model* model){ + + extern int my_rank; + + /*fem model: */ + FemModel* fem_model=NULL; + + /*output: */ + DataSet* results=NULL; + Result* result=NULL; + + /*Intermediary: */ + DataSet* diagnostic_results=NULL; + DataSet* gradjcompute_results=NULL; + DataSet* steadystate_results=NULL; + Vec u_g=NULL; + Vec t_g=NULL; + Vec m_g=NULL; + double search_scalar; + int control_type; + double* fit=NULL; + double* optscal=NULL; + int gsize; + double* maxiter=NULL; + double* cm_jump=NULL; + double eps_cm; + double tolx; + double* param_g=NULL; + Vec grad_g=NULL; + Vec new_grad_g=NULL; + Vec grad_g_old=NULL; + double* grad_g_double=NULL; + double cm_min; + double cm_max; + int cm_gradient; + int nsteps,n,i; + double* J=NULL; + OptArgs optargs; + OptPars optpars; + Param* param=NULL; + + /*flags: */ + int analysis_type; + int sub_analysis_type; + int control_steady; + int verbose=0; + int converged=0; + int numberofnodes; + + //initialize results + results=new DataSet(ResultsEnum); + + /*Process models*/ + ControlInitialization(model); + fem_model=model->GetActiveFormulation(); + + /*Recover parameters used throughout the solution:*/ + model->FindParam(&nsteps,NStepsEnum); + model->FindParam(&control_type,ControlTypeEnum); + model->FindParam(&fit,NULL,NULL,FitEnum); + model->FindParam(&optscal,NULL,NULL,OptScalEnum); + model->FindParam(&maxiter,NULL,NULL,MaxIterEnum); + model->FindParam(&cm_jump,NULL,NULL,CmJumpEnum); + model->FindParam(&eps_cm,EpsCmEnum); + model->FindParam(&tolx,TolXEnum); + model->FindParam(&cm_min,CmMinEnum); + model->FindParam(&cm_max,CmMaxEnum); + model->FindParam(&cm_gradient,CmGradientEnum); + model->FindParam(¶m_g,NULL,NULL,ControlParameterEnum); + model->FindParam(&analysis_type,AnalysisTypeEnum); + model->FindParam(&sub_analysis_type,SubAnalysisTypeEnum); + model->FindParam(&numberofnodes,NumberOfNodesEnum); + model->FindParam(&control_steady,ControlSteadyEnum); + gsize=fem_model->nodes->NumberOfDofs(); + + /*Initialize misfit: */ + J=(double*)xmalloc(nsteps*sizeof(double)); + + /*Start looping: */ + for(n=0;nUpdateInputsFromVector(param_g,control_type,VertexEnum); + model->UpdateInputsFromConstant(fit[n],FitEnum); + + /*In case we are running a steady state control method, compute new temperature field using new parameter + * distribution: */ + if (control_steady){ + steadystate_results= steadystate_core(model); + VecFree(&t_g); steadystate_results->FindResult(&t_g,"t_g"); + delete steadystate_results; + model->UpdateInputsFromVector(t_g,TemperatureEnum,VertexEnum); + } + + _printf_("%s\n"," computing gradJ..."); + gradjcompute_results= gradjcompute_core(model); + gradjcompute_results->FindResult(&grad_g,"grad_g"); + delete gradjcompute_results; + + /*Return gradient if asked: */ + if (cm_gradient){ + + /*Plug results into output dataset: */ + result=new Result(results->Size()+1,0,1,"grad_g",grad_g); + results->AddObject(result); + + result=new Result(results->Size()+1,0,1,"analysis_type",EnumAsString(analysis_type)); + results->AddObject(result); + + /*Free ressources: */ + xfree((void**)&control_type); + xfree((void**)&fit); + xfree((void**)&optscal); + xfree((void**)&maxiter); + xfree((void**)&cm_jump); + xfree((void**)&grad_g_double); + xfree((void**)¶m_g); + VecFree(&u_g); + VecFree(&t_g); + VecFree(&m_g); + xfree((void**)&J); + return results; + } + + /*Normalize if last gradient not satisfying (search_scalar==0)*/ + if (n>0 && search_scalar==0){ + _printf_("%s"," orthogonalization..."); + Orthx(&new_grad_g,grad_g,grad_g_old); + _printf_("%s\n"," done."); + } + else{ + _printf_("%s"," normalizing directions..."); + Orthx(&new_grad_g,grad_g,NULL); + _printf_("%s\n"," done."); + } + VecFree(&grad_g); VecFree(&grad_g_old); + grad_g_old=new_grad_g; + VecToMPISerial(&grad_g_double,new_grad_g); + + _printf_("%s\n"," optimizing along gradient direction"); + optargs.model=model; + optargs.param_g=param_g; optargs.grad_g=grad_g_double; optargs.n=n; + optpars.xmin=0; optpars.xmax=1; optpars.tolerance=tolx; optpars.maxiter=(int)maxiter[n];optpars.cm_jump=cm_jump[n]; + BrentSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs); + + _printf_("%s"," updating parameter using optimized search scalar..."); + for(i=0;i=0){ + if (fit[i]==fit[n]){ + //convergence test only if we have the same misfits + if ((J[i]-J[n])/J[n] <= eps_cm){ + //convergence if convergence criteria fullfilled + converged=1; + _printf_("%s%g%s%g\n"," Convergence criterion: dJ/J = ",(J[i]-J[n])/J[n],"<",eps_cm); + } + else{ + _printf_("%s%g%s%g\n"," Convergence criterion: dJ/J = ",(J[i]-J[n])/J[n],">",eps_cm); + } + break; + } + i=i-1; + } + } + //stop if convergence has been reached + if(converged) break; + + //some temporary saving + if (((n+1)%5)==0){ + _printf_("%s"," saving temporary results..."); + ControlRestart(model,param_g); + _printf_("%s\n"," done."); + } + } + + /*Write results to disk: */ + _printf_("%s"," preparing final velocity solution"); + /*Launch diagnostic with the last parameter distribution*/ + if (control_steady){ + model->UpdateInputsFromVector(param_g,control_type,VertexEnum); + steadystate_results=steadystate_core(model); + + //extract u_g ,t_g and m_g from steadystate results, and erase diagnostic_results; + steadystate_results->FindResult(&u_g,"u_g"); + steadystate_results->FindResult(&m_g,"m_g"); + steadystate_results->FindResult(&t_g,"t_g"); + delete steadystate_results; + } + else{ + model->UpdateInputsFromVector(param_g,control_type,VertexEnum); + diagnostic_results=diagnostic_core(model); + + //extract u_g from diagnostic_results, and erase diagnostic_results; + diagnostic_results->FindResult(&u_g,"u_g"); + delete diagnostic_results; + } + + /*Plug results into output dataset: */ + result=new Result(results->Size()+1,0,1,"u_g",u_g); + results->AddObject(result); + result=new Result(results->Size()+1,0,1,"param_g",param_g,numberofnodes); + results->AddObject(result); + result=new Result(results->Size()+1,0,1,"J",J,nsteps); + results->AddObject(result); + if (control_steady){ + result=new Result(results->Size()+1,0,1,"t_g",t_g); + results->AddObject(result); + result=new Result(results->Size()+1,0,1,"m_g",m_g); + results->AddObject(result); + } + + /*Add analysis_type and control_type to results: */ + result=new Result(results->Size()+1,0,1,"analysis_type",EnumAsString(analysis_type)); + results->AddObject(result); + result=new Result(results->Size()+1,0,1,"control_type",EnumAsString(control_type)); + results->AddObject(result); + + + /*Free ressources: */ + xfree((void**)&control_type); + xfree((void**)&fit); + xfree((void**)&optscal); + xfree((void**)&maxiter); + xfree((void**)&cm_jump); + VecFree(&new_grad_g); //do not VecFree grad_g and grad_g_old, they point to new_grad_g + xfree((void**)&grad_g_double); + xfree((void**)¶m_g); + VecFree(&u_g); + VecFree(&t_g); + VecFree(&m_g); + xfree((void**)&J); + + //return: + return results; +} Index: /issm/trunk/src/c/solutions/convergence.cpp =================================================================== --- /issm/trunk/src/c/solutions/convergence.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/convergence.cpp (revision 3894) @@ -0,0 +1,141 @@ +/*!\file: convergence.cpp + * \brief: figure out if convergence has been reached + */ + +#include "../objects/objects.h" +#include "../modules.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +void convergence(int* pconverged, Mat Kff,Vec pf,Vec uf,Vec old_uf,Parameters* parameters){ + + /*output*/ + int converged=0; + + /*intermediary*/ + Vec KU=NULL; + Vec KUF=NULL; + Vec KUold=NULL; + Vec KUoldF=NULL; + Vec duf=NULL; + double ndu,nduinf,nu; + double nKUF; + double nKUoldF; + double nF; + double solver_residue,res; + + /*convergence options*/ + double eps_res; + double eps_rel; + double eps_abs; + int verbose; + double yts; + + /*Check arguments*/ + ISSMASSERT(Kff && pf && uf && old_uf); + + /*get convergence options*/ + parameters->FindParam(&eps_res,EpsResEnum); + parameters->FindParam(&eps_rel,EpsRelEnum); + parameters->FindParam(&eps_abs,EpsAbsEnum); + parameters->FindParam(&yts,YtsEnum); + parameters->FindParam(&verbose,VerboseEnum); + + /*Display solver caracteristics*/ + if (verbose>1){ + + //compute KUF = KU - F = K*U - F + VecDuplicate(uf,&KU); MatMultPatch(Kff,uf,KU); + VecDuplicate(KU,&KUF);VecCopy(KU,KUF); VecAYPX(KUF,-1.0,pf); + + //compute norm(KUF), norm(F) and residue + VecNorm(KUF,NORM_2,&nKUF); + VecNorm(pf,NORM_2,&nF); + solver_residue=nKUF/nF; + _printf_("\n%s%g\n"," solver residue: norm(KU-F)/norm(F)=",solver_residue); + + //clean up + VecFree(&KU); + VecFree(&KUF); + } + + /*Force equilibrium (Mandatory)*/ + + //compute K[n]U[n-1]F = K[n]U[n-1] - F + VecDuplicate(uf,&KUold); MatMultPatch(Kff,old_uf,KUold); + VecDuplicate(KUold,&KUoldF);VecCopy(KUold,KUoldF); VecAYPX(KUoldF,-1.0,pf); + VecNorm(KUoldF,NORM_2,&nKUoldF); + VecNorm(pf,NORM_2,&nF); + res=nKUoldF/nF; + if (isnan(res)){ + _printf_("norm nf = %lf and norm kuold = %lf\n",nF,nKUoldF); + ISSMERROR("mechanical equilibrium convergence criterion is NaN!"); + } + + //clean up + VecFree(&KUold); + VecFree(&KUoldF); + + //print + if(res ",eps_res*100," \%"); + converged=0; + } + + /*Relative criterion (optional)*/ + if (!isnan(eps_rel) || (verbose>1)){ + + //compute norm(du)/norm(u) + VecDuplicate(old_uf,&duf);VecCopy(old_uf,duf); VecAYPX(duf,-1.0,uf); + VecNorm(duf,NORM_2,&ndu); VecNorm(old_uf,NORM_2,&nu); + + if (isnan(ndu) || isnan(nu)) ISSMERROR("convergence criterion is NaN!"); + + //clean up + VecFree(&duf); + + //print + if (!isnan(eps_rel)){ + if((ndu/nu) ",eps_rel*100," \%"); + converged=0; + } + } + else _printf_("%-50s%g%s\n"," Convergence criterion: norm(du)/norm(u)",ndu/nu*100," \%"); + + } + + /*Absolute criterion (Optional) = max(du)*/ + if (!isnan(eps_abs) || (verbose>1)){ + + //compute max(du) + VecDuplicate(old_uf,&duf);VecCopy(old_uf,duf); VecAYPX(duf,-1.0,uf); + VecNorm(duf,NORM_2,&ndu); VecNorm(duf,NORM_INFINITY,&nduinf); + if (isnan(ndu) || isnan(nu)) ISSMERROR("convergence criterion is NaN!"); + + //clean up + VecFree(&duf); + + //print + if (!isnan(eps_abs)){ + if ((nduinf*yts) ",eps_abs," m/yr"); + converged=0; + } + } + else _printf_("%-50s%g%s\n"," Convergence criterion: max(du)",nduinf*yts," m/yr"); + + } + + /*assign output*/ + *pconverged=converged; +} Index: /issm/trunk/src/c/solutions/diagnostic.cpp =================================================================== --- /issm/trunk/src/c/solutions/diagnostic.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/diagnostic.cpp (revision 3894) @@ -0,0 +1,154 @@ +/*!\file: diagnostic.cpp + * \brief: diagnostic solution + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../DataSet/DataSet.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../include/include.h" +#include "../modules.h" +#include "./parallel.h" + +int main(int argc,char* *argv){ + + /*I/O: */ + FILE* fid=NULL; + char* inputfilename=NULL; + char* outputfilename=NULL; + char* lockname=NULL; + bool qmu_analysis=false; + bool control_analysis=false; + + /*Model: */ + Model* model=NULL; + + /*Results: */ + DataSet* results=NULL; + DataSet* processed_results=NULL; + Result* result=NULL; + + bool waitonlock=false; + + /*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"); + + _printf_("read and create finite element model:\n"); + _printf_("\n reading diagnostic horiz model data:\n"); + model->AddFormulation(fid,DiagnosticAnalysisEnum,HorizAnalysisEnum); + + _printf_("\n reading diagnostic vert model data:\n"); + model->AddFormulation(fid,DiagnosticAnalysisEnum,VertAnalysisEnum); + + _printf_("\n reading diagnostic stokes model data:\n"); + model->AddFormulation(fid,DiagnosticAnalysisEnum,StokesAnalysisEnum); + + _printf_("\n reading diagnostic hutter model data:\n"); + model->AddFormulation(fid,DiagnosticAnalysisEnum,HutterAnalysisEnum); + + _printf_("\n reading surface and bed slope computation model data:\n"); + model->AddFormulation(fid,SlopecomputeAnalysisEnum); + + /*get parameters: */ + model->FindParam(&qmu_analysis,QmuAnalysisEnum); + model->FindParam(&control_analysis,ControlAnalysisEnum); + model->FindParam(&waitonlock,WaitOnLockEnum); + + MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime(); + + /*are we running the solution sequence, or a qmu wrapper around it? : */ + if(!qmu_analysis){ + if(!control_analysis){ + + _printf_("call computational core:\n"); + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + results=diagnostic_core(model); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + + _printf_("process results:\n"); + ProcessResults(&processed_results,results,model,DiagnosticAnalysisEnum); + } + else{ + /*run control analysis: */ + _printf_("call computational core:\n"); + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + results=control_core(model); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + + _printf_("process results:\n"); + ProcessResults(&processed_results,results,model,ControlAnalysisEnum); + } + + _printf_("write results to disk:\n"); + OutputResults(processed_results,outputfilename); + } + else{ + /*run qmu analysis: */ + _printf_("calling qmu analysis on diagnostic core:\n"); + + #ifdef _HAVE_DAKOTA_ + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + Qmux(model,DiagnosticAnalysisEnum,NoneAnalysisEnum); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + #else + ISSMERROR(" Dakota not present, cannot do qmu!"); + #endif + } + + if (waitonlock>0){ + _printf_("write lock file:\n"); + WriteLockFile(lockname); + } + + /*Free ressources */ + delete model; + delete results; + delete processed_results; + + /*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/solutions/diagnostic_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/diagnostic_core.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/diagnostic_core.cpp (revision 3894) @@ -0,0 +1,209 @@ +/*!\file: diagnostic_core.cpp + * \brief: core of the diagnostic solution + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./parallel.h" +#include "../modules.h" +#include "../include/include.h" + +DataSet* diagnostic_core(Model* model){ + + extern int my_rank; + int dummy; + + /*fem models: */ + FemModel* fem_dh=NULL; + FemModel* fem_dv=NULL; + FemModel* fem_dhu=NULL; + FemModel* fem_ds=NULL; + FemModel* fem_sl=NULL; + + /*output: */ + DataSet* results=NULL; + Result* result=NULL; + + /*solutions: */ + Vec ug=NULL; + Vec ug_horiz=NULL; + Vec ug_vert=NULL; + Vec ug_stokes=NULL; + Vec pg=NULL; + Vec riftproperties=NULL; + double* u_g_initial=NULL; + double* vx=NULL; + double* vy=NULL; + double* vz=NULL; + + /*flags: */ + int verbose=0; + bool qmu_analysis=false; + int dim=-1; + bool ishutter=false; + bool ismacayealpattyn=false; + bool isstokes=false; + int numberofdofspernode_sl; + int numberofdofspernode_dh; + int numberofdofspernode_ds; + int numberofnodes; + int numrifts=0; + + /*slopes: */ + Vec surfaceslopex=NULL; + Vec surfaceslopey=NULL; + Vec bedslopex=NULL; + Vec bedslopey=NULL; + + double stokesreconditioning; + + /*dof recovery: */ + int dof01[2]={0,1}; + int dof2[1]={2}; + int dof012[3]={0,1,2}; + int dof3[1]={3}; + double* dofset=NULL; + + //initialize results + results=new DataSet(ResultsEnum); + + //first recover parameters common to all solutions + model->FindParam(&verbose,VerboseEnum); + model->FindParam(&dim,DimEnum); + model->FindParam(&ishutter,IsHutterEnum); + model->FindParam(&ismacayealpattyn,IsMacAyealPattynEnum); + model->FindParam(&numberofnodes,NumberOfNodesEnum); + model->FindParam(&isstokes,IsStokesEnum); + model->FindParam(&stokesreconditioning,StokesReconditioningEnum); + model->FindParam(&numrifts,NumRiftsEnum); + model->FindParam(&qmu_analysis,QmuAnalysisEnum); + + /*recover fem models: */ + fem_dh=model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum); + fem_dv=model->GetFormulation(DiagnosticAnalysisEnum,VertAnalysisEnum); + fem_ds=model->GetFormulation(DiagnosticAnalysisEnum,StokesAnalysisEnum); + fem_dhu=model->GetFormulation(DiagnosticAnalysisEnum,HutterAnalysisEnum); + fem_sl=model->GetFormulation(SlopecomputeAnalysisEnum); + + //specific parameters for specific models + fem_dh->FindParam(&numberofdofspernode_dh,NumberOfDofsPerNodeEnum); + fem_sl->FindParam(&numberofdofspernode_sl,NumberOfDofsPerNodeEnum); + fem_ds->FindParam(&numberofdofspernode_ds,NumberOfDofsPerNodeEnum); + + //for qmu analysis, be sure the velocity input we are starting from is the one in the parameters: */ + if(qmu_analysis){ + fem_dh->FindParam(&vx,&dummy,VxEnum); model->UpdateInputsFromVector(vx,VxEnum,VertexEnum); + fem_dh->FindParam(&vy,&dummy,VyEnum); model->UpdateInputsFromVector(vy,VyEnum,VertexEnum); + fem_dh->FindParam(&vz,&dummy,VzEnum); model->UpdateInputsFromVector(vz,VzEnum,VertexEnum); + } + + /*Compute slopes: */ + slope_core(&surfaceslopex,&surfaceslopey,fem_sl,SurfaceAnalysisEnum); + slope_core(&bedslopex,&bedslopey,fem_sl,BedAnalysisEnum); + + /*Update: */ + model->UpdateInputsFromVector(surfaceslopex,SurfaceSlopexEnum,VertexEnum); + model->UpdateInputsFromVector(surfaceslopey,SurfaceSlopeyEnum,VertexEnum); + model->UpdateInputsFromVector(bedslopex,BedSlopexEnum,VertexEnum); + model->UpdateInputsFromVector(bedslopey,BedSlopeyEnum,VertexEnum); + + + if(ishutter){ + + if(verbose)_printf_("%s\n"," computing hutter velocities..."); + diagnostic_core_linear(&ug,fem_dhu,DiagnosticAnalysisEnum,HutterAnalysisEnum); + + if(verbose)_printf_("%s\n"," computing pressure according to MacAyeal..."); + ComputePressurex(&pg,fem_dhu->elements,fem_dhu->nodes, fem_dhu->vertices,fem_dhu->loads,fem_dhu->materials,fem_dhu->parameters,DiagnosticAnalysisEnum,HutterAnalysisEnum); + + if(verbose)_printf_("%s\n"," update boundary conditions for macyeal pattyn using hutter results..."); + if (ismacayealpattyn){ + VecFree(&fem_dh->yg->vector); VecFree(&fem_dh->ys);VecFree(&fem_dh->ys0); + VecDuplicatePatch(&fem_dh->yg->vector,ug); + Reducevectorgtosx(&fem_dh->ys,&fem_dh->ys0, fem_dh->yg->vector,fem_dh->nodesets); + } + + } + + if (ismacayealpattyn){ + + if(verbose)_printf_("%s\n"," computing horizontal velocities..."); + diagnostic_core_nonlinear(&ug,NULL,NULL,fem_dh->loads,fem_dh,DiagnosticAnalysisEnum,HorizAnalysisEnum); + + if(dim==2){ + if(verbose)_printf_("%s\n"," computing pressure according to MacAyeal..."); + ComputePressurex(&pg,fem_dh->elements,fem_dh->nodes, fem_dh->vertices,fem_dh->loads,fem_dh->materials,fem_dh->parameters,DiagnosticAnalysisEnum,HorizAnalysisEnum); + } + + } + + + if (dim==3){ + + if(verbose)_printf_("%s\n"," extruding horizontal velocities..."); + FieldExtrudex(ug,fem_dh->elements,fem_dh->nodes, fem_dh->vertices,fem_dh->loads,fem_dh->materials,fem_dh->parameters,"velocity",1); + SplitSolutionVectorx(ug,numberofnodes,numberofdofspernode_dh,&vx,&vy); + model->UpdateInputsFromVector(vx,VxEnum,VertexEnum); + model->UpdateInputsFromVector(vy,VyEnum,VertexEnum); + + if(verbose)_printf_("%s\n"," computing vertical velocities..."); + diagnostic_core_linear(&ug_vert,fem_dv,DiagnosticAnalysisEnum,VertAnalysisEnum); + model->UpdateInputsFromVector(ug_vert,VzEnum,VertexEnum); + + if(verbose)_printf_("%s\n"," computing pressure according to Pattyn..."); + ComputePressurex(&pg,fem_dh->elements, fem_dh->nodes, fem_dh->vertices,fem_dh->loads, fem_dh->materials,fem_dh->parameters,DiagnosticAnalysisEnum,HorizAnalysisEnum); + model->UpdateInputsFromVector(pg,PressureEnum,VertexEnum); + + if (isstokes){ + + //"recondition" pressure + VecScale(pg,1.0/stokesreconditioning); + + if(verbose)_printf_("%s\n"," update boundary conditions for stokes using velocities previously computed..."); + xfree((void**)&dofset);dofset=dofsetgen(2,dof01,4,numberofnodes*4); VecMerge(fem_ds->yg->vector,ug,dofset,2*numberofnodes); + xfree((void**)&dofset);dofset=dofsetgen(1,dof2,4,numberofnodes*4); VecMerge(fem_ds->yg->vector,ug_vert,dofset,1*numberofnodes); + VecFree(&fem_ds->ys); VecFree(&fem_ds->ys0); + Reducevectorgtosx(&fem_ds->ys,&fem_ds->ys0, fem_ds->yg->vector,fem_ds->nodesets); + + if(verbose)_printf_("%s\n"," computing stokes velocities and pressure ..."); + VecFree(&ug); + diagnostic_core_nonlinear(&ug,NULL,NULL,NULL,fem_ds,DiagnosticAnalysisEnum,StokesAnalysisEnum); + + //decondition" pressure + VecFree(&pg); + xfree((void**)&dofset);dofset=dofsetgen(1,dof3,4,numberofnodes*4); VecPartition(&pg, ug, dofset, numberofnodes*1); + VecScale(pg,stokesreconditioning); + } + } + + /*Plug results into output dataset: */ + result=new Result(results->Size()+1,0,1,"u_g",ug); + results->AddObject(result); + result=new Result(results->Size()+1,0,1,"p_g",pg); + results->AddObject(result); + result=new Result(results->Size()+1,0,1,"analysis_type",EnumAsString(DiagnosticAnalysisEnum)); + results->AddObject(result); + + /*output if we have rifts: */ + if(numrifts){ + OutputRiftsx( &riftproperties,fem_dh->loads,numrifts); + result=new Result(results->Size()+1,0,1,"riftproperties",riftproperties); + results->AddObject(result); + } + + /*Free ressources: */ + VecFree(&ug); + VecFree(&pg); + xfree((void**)&dofset); + xfree((void**)&u_g_initial); + + VecFree(&surfaceslopex); + VecFree(&surfaceslopey); + VecFree(&bedslopex); + VecFree(&bedslopey); + + /*return: */ + return results; +} Index: /issm/trunk/src/c/solutions/diagnostic_core_linear.cpp =================================================================== --- /issm/trunk/src/c/solutions/diagnostic_core_linear.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/diagnostic_core_linear.cpp (revision 3894) @@ -0,0 +1,62 @@ +/*!\file: diagnostic_core_nonlinear.cpp + * \brief: core of the diagnostic solution for non linear materials + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../modules.h" + +void diagnostic_core_linear(Vec* pug,FemModel* fem,int analysis_type,int sub_analysis_type){ + + /*parameters:*/ + int kflag,pflag; + int verbose=0; + char* solver_string=NULL; + + /*output: */ + Vec ug=NULL; + Vec uf=NULL; + + /*intermediary: */ + Mat Kgg=NULL; + Mat Kff=NULL; + Mat Kfs=NULL; + Vec pg=NULL; + Vec pf=NULL; + + /*Recover parameters: */ + kflag=1; pflag=1; + fem->parameters->FindParam(&verbose,VerboseEnum); + fem->parameters->FindParam(&solver_string,SolverStringEnum); + + //*Generate system matrices + if (verbose) _printf_(" Generating matrices\n"); + SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); + + if (verbose) _printf_(" Generating penalty matrices\n"); + //*Generate penalty system matrices + PenaltySystemMatricesx(Kgg, pg,NULL,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); + + /*!Reduce matrix from g to f size:*/ + if (verbose) _printf_(" reducing matrix from g to f set\n"); + Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets); MatFree(&Kgg); + + /*!Reduce load from g to f size: */ + if (verbose) _printf_(" reducing load from g to f set\n"); + Reduceloadfromgtofx(&pf, pg, fem->Gmn, Kfs, fem->ys, fem->nodesets);VecFree(&pg); MatFree(&Kfs); + + /*Solve: */ + if (verbose) _printf_(" solving\n"); + Solverx(&uf, Kff, pf, NULL, solver_string); MatFree(&Kff); VecFree(&pf); + + //Merge back to g set + if (verbose) _printf_(" merging solution from f to g set\n"); + Mergesolutionfromftogx(&ug, uf,fem->Gmn,fem->ys,fem->nodesets);VecFree(&uf); + + /*free ressources: */ + xfree((void**)&solver_string); + + /*Assign output pointers:*/ + *pug=ug; +} Index: /issm/trunk/src/c/solutions/diagnostic_core_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/solutions/diagnostic_core_nonlinear.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/diagnostic_core_nonlinear.cpp (revision 3894) @@ -0,0 +1,158 @@ +/*!\file: diagnostic_core_nonlinear.cpp + * \brief: core of the diagnostic solution for non linear materials + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../modules.h" +#include "./parallel.h" + +void diagnostic_core_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, DataSet* input_loads,FemModel* fem,int analysis_type,int sub_analysis_type){ + + + /*solution : */ + Vec ug=NULL; + Vec uf=NULL; + Vec old_ug=NULL; + Vec old_uf=NULL; + DataSet* loads=NULL; + + /*intermediary: */ + Mat Kgg=NULL; + Mat Kff=NULL; + Mat Kfs=NULL; + Vec pg=NULL; + Vec pf=NULL; + int converged; + int constraints_converged; + int num_unstable_constraints; + int count; + int numberofnodes; + int min_mechanical_constraints; + int max_nonlinear_iterations; + + /*parameters:*/ + int kflag,pflag; + char* solver_string=NULL; + int verbose=0; + + /*Recover parameters: */ + kflag=1; pflag=1; + fem->FindParam(&numberofnodes,NumberOfNodesEnum); + fem->FindParam(&solver_string,SolverStringEnum); + fem->FindParam(&verbose,VerboseEnum); + fem->FindParam(&min_mechanical_constraints,MinMechanicalConstraintsEnum); + fem->FindParam(&max_nonlinear_iterations,MaxNonlinearIterationsEnum); + + /*Were loads requested as output? : */ + if(!input_loads){ + loads=fem->loads->Copy(); //we don't want to clobber loads, as they are not needed in output. + } + else{ + loads=input_loads; //we are going to modify the loads! + } + + count=1; + converged=0; + + /*Start non-linear iteration using input velocity: */ + GetSolutionFromInputsx(&ug, fem->elements, fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters, analysis_type, sub_analysis_type); + Reducevectorgtofx(&uf, ug, fem->nodesets); + + for(;;){ + + //save pointer to old velocity + VecFree(&old_ug);old_ug=ug; + VecFree(&old_uf);old_uf=uf; + + if (verbose) _printf_(" Generating matrices\n"); + //*Generate system matrices + SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); + + if (verbose) _printf_(" Generating penalty matrices\n"); + //*Generate penalty system matrices + PenaltySystemMatricesx(Kgg, pg,NULL,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); + + if (verbose) _printf_(" reducing matrix from g to f set\n"); + /*!Reduce matrix from g to f size:*/ + Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets); + + /*Free ressources: */ + MatFree(&Kgg); + + if (verbose) _printf_(" reducing load from g to f set\n"); + /*!Reduce load from g to f size: */ + Reduceloadfromgtofx(&pf, pg, fem->Gmn, Kfs, fem->ys, fem->nodesets); + + //no need for pg and Kfs anymore + VecFree(&pg); + MatFree(&Kfs); + + /*Solve: */ + if (verbose) _printf_(" solving\n"); + Solverx(&uf, Kff, pf, old_uf, solver_string); + + //Merge back to g set + if (verbose) _printf_(" merging solution from f to g set\n"); + Mergesolutionfromftogx(&ug, uf,fem->Gmn,fem->ys,fem->nodesets); + + //Update inputs using new solution: + fem->UpdateInputsFromSolution(ug,analysis_type,sub_analysis_type); + + //Deal with penalty loads + if (verbose) _printf_(" penalty constraints\n"); + PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,analysis_type,sub_analysis_type); + + if(verbose)_printf_(" number of unstable constraints: %i\n",num_unstable_constraints); + + /*Figure out if convergence is reached.*/ + convergence(&converged,Kff,pf,uf,old_uf,fem->parameters); + MatFree(&Kff);VecFree(&pf); + + /*add converged to inputs: */ + fem->UpdateInputsFromConstant(converged,ConvergedEnum); + + //rift convergence + if (!constraints_converged) { + if (converged){ + if (num_unstable_constraints <= min_mechanical_constraints) converged=1; + else converged=0; + } + } + + /*Increase count: */ + count++; + if(converged==1)break; + if(count>=max_nonlinear_iterations){ + _printf_(" maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations); + break; + } + } + + //more output might be needed, when running in control + if(pKff0){ + + kflag=1; pflag=0; //stiffness generation only + + SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); + Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets); + MatFree(&Kgg);VecFree(&pg); + + } + + /*Delete loads only if no ouput was requested: */ + if(!input_loads)delete loads; + + /*clean up*/ + VecFree(&uf); + VecFree(&old_uf); + VecFree(&old_ug); + xfree((void**)&solver_string); + + /*Assign output pointers: */ + *pug=ug; + if(pKff0)*pKff0=Kff; + if(pKfs0)*pKfs0=Kfs; + +} Index: /issm/trunk/src/c/solutions/gradjcompute_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/gradjcompute_core.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/gradjcompute_core.cpp (revision 3894) @@ -0,0 +1,135 @@ +/*!\file: gradjcompute_core.cpp + * \brief compute inverse method gradient direction. + */ + +#include "../modules.h" +#include "./parallel.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +DataSet* gradjcompute_core(Model* model){ + + + /*intermediary: */ + FemModel* femmodel=NULL; + DataSet* diagnostic_results=NULL; + int analysis_type; + int sub_analysis_type; + int numberofnodes; + int control_steady; + int numberofdofspernode; + char* solverstring=NULL; + int control_type; + + Vec u_g=NULL; + + Vec du_g=NULL; + Vec du_f=NULL; + Vec grad_g=NULL; + + Vec lambda_f=NULL; + Vec lambda_g=NULL; + double* lambdax=NULL; + double* lambday=NULL; + double* vx=NULL; + double* vy=NULL; + double* vz=NULL; + + Mat K_ff0=NULL; + Mat K_fs0=NULL; + + /*output: */ + DataSet* results=NULL; + Result* result=NULL; + + /*flags: */ + int verbose=0; + int dim=-1; + int extrude_param=0; + + //initialize results + results=new DataSet(ResultsEnum); + + /*some parameters:*/ + femmodel=model->GetActiveFormulation(); + femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum); + femmodel->parameters->FindParam(&sub_analysis_type,SubAnalysisTypeEnum); + femmodel->parameters->FindParam(&control_steady,ControlSteadyEnum); + femmodel->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + femmodel->parameters->FindParam(&numberofdofspernode,NumberOfDofsPerNodeEnum); + femmodel->parameters->FindParam(&solverstring,SolverStringEnum); + femmodel->parameters->FindParam(&control_type,ControlTypeEnum); + femmodel->parameters->FindParam(&extrude_param,ExtrudeParamEnum); + femmodel->parameters->FindParam(&verbose,VerboseEnum); + femmodel->parameters->FindParam(&dim,DimEnum); + + _printf_("%s\n"," recover solution for this stiffness and right hand side:"); + diagnostic_core_nonlinear(&u_g,&K_ff0,&K_fs0,NULL, femmodel,DiagnosticAnalysisEnum,sub_analysis_type); + + _printf_("%s\n"," buid Du, difference between observed velocity and model velocity:"); + Dux( &du_g, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,analysis_type,sub_analysis_type); + + _printf_("%s\n"," reduce adjoint load from g-set to f-set:"); + Reduceloadfromgtofx(&du_f, du_g, femmodel->Gmn, K_fs0, femmodel->ys0, femmodel->nodesets); + VecFree(&du_g);MatFree(&K_fs0); + + _printf_("%s\n"," solve for adjoint vector:"); + Solverx(&lambda_f, K_ff0, du_f, NULL, solverstring); + VecFree(&du_f); MatFree(&K_ff0); + + _printf_("%s\n"," merge back to g set:"); + Mergesolutionfromftogx(&lambda_g, lambda_f,femmodel->Gmn,femmodel->ys0,femmodel->nodesets); + VecFree(&lambda_f); + + /*add to inputs: */ + SplitSolutionVectorx(lambda_g,numberofnodes,numberofdofspernode,&lambdax,&lambday); + model->UpdateInputsFromVector(lambdax,AdjointxEnum,VertexEnum); + model->UpdateInputsFromVector(lambday,AdjointyEnum,VertexEnum); + VecFree(&lambda_g); + + _printf_("%s\n"," compute gradJ:"); + Gradjx( &grad_g, numberofnodes,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, + analysis_type,sub_analysis_type,control_type); + + if (dim==3 && extrude_param){ + + _printf_("%s\n"," extruding gradient..."); + FieldExtrudex( grad_g, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"gradj",0); + } + + if(control_steady){ + diagnostic_results= diagnostic_core(model); + + //extract u_g and add it to input (3d velocity needed by thermal_core) + diagnostic_results->FindResult(&u_g,"u_g"); + + SplitSolutionVectorx(u_g,numberofnodes,3,&vx,&vy,&vz); + model->UpdateInputsFromVector(vx,VxEnum,VertexEnum); + model->UpdateInputsFromVector(vy,VyEnum,VertexEnum); + model->UpdateInputsFromVector(vz,VzEnum,VertexEnum); + + delete diagnostic_results; + } + + /*Plug results into output dataset: */ + result=new Result(results->Size()+1,0,1,"grad_g",grad_g); + results->AddObject(result); + result=new Result(results->Size()+1,0,1,"analysis_type",EnumAsString(GradientAnalysisEnum)); + results->AddObject(result); + + /*Free ressources:*/ + VecFree(&u_g); + VecFree(&grad_g); + xfree((void**)&solverstring); + xfree((void**)&lambdax); + xfree((void**)&lambday); + xfree((void**)&vx); + xfree((void**)&vy); + xfree((void**)&vz); + +} Index: /issm/trunk/src/c/solutions/objectivefunctionC.cpp =================================================================== --- /issm/trunk/src/c/solutions/objectivefunctionC.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/objectivefunctionC.cpp (revision 3894) @@ -0,0 +1,127 @@ +/*!\file: objectivefunctionC + * \brief objective function that returns a misfit, for a certain parameter. + */ + +#include "../modules.h" +#include "./parallel.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +double objectivefunctionC(double search_scalar,OptArgs* optargs){ + + int i; + + /*output: */ + double J; + + /*parameters: */ + Model* model=NULL; + FemModel* femmodel=NULL; + DataSet* diagnostic_results=NULL; + double* param_g=NULL; + double* grad_g=NULL; + int numberofdofspernode; + int n; + + /*intermediary:*/ + int gsize; + double* optscal=NULL; + double* fit=NULL; + double cm_min; + double cm_max; + int control_type; + double* param_g_copy=NULL; + int analysis_type; + int sub_analysis_type; + int control_steady; + Vec u_g=NULL; + Vec u_g_full=NULL; + double* u_g_double=NULL; + double* vx=NULL; + double* vy=NULL; + double* vz=NULL; + int numberofnodes; + + /*steadystate: */ + int dt=0; + int isstokes=0; + DataSet* results_steadystate=NULL; + int dofs01[2]={0,1}; + double* dofset=NULL; + + /*Recover active model: */ + model=optargs->model; + femmodel=model->GetActiveFormulation(); + + /*Recover parameters: */ + param_g=optargs->param_g; + grad_g=optargs->grad_g; + n=optargs->n; + + gsize=femmodel->nodesets->GetGSize(); + femmodel->parameters->FindParam(&optscal,NULL,NULL,OptScalEnum); + femmodel->parameters->FindParam(&fit,NULL,NULL,FitEnum); + femmodel->parameters->FindParam(&cm_min,CmMinEnum); + femmodel->parameters->FindParam(&cm_max,CmMaxEnum); + femmodel->parameters->FindParam(&control_type,ControlTypeEnum); + femmodel->parameters->FindParam(&control_steady,ControlSteadyEnum); + femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum); + femmodel->parameters->FindParam(&sub_analysis_type,SubAnalysisTypeEnum); + femmodel->parameters->FindParam(&isstokes,IsStokesEnum); + femmodel->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + femmodel->parameters->FindParam(&numberofdofspernode,NumberOfDofsPerNodeEnum); + + /*First copy param_g so we don't modify it: */ + param_g_copy=(double*)xmalloc(numberofnodes*sizeof(double)); + memcpy(param_g_copy,param_g,numberofnodes*sizeof(double)); + + /*First, update param_g using search_scalar: */ + for(i=0;iUpdateInputsFromVector(param_g_copy,control_type,VertexEnum); + + /*Run diagnostic with updated parameters.*/ + if(!control_steady){ + diagnostic_core_nonlinear(&u_g,NULL,NULL,NULL,femmodel,DiagnosticAnalysisEnum,sub_analysis_type); + femmodel->UpdateInputsFromSolution(u_g,DiagnosticAnalysisEnum,sub_analysis_type); + VecFree(&u_g); + } + else{ + //We need a 3D velocity!! (vz is required for the next thermal run) + diagnostic_results= diagnostic_core(model); + + //extract u_g and add it to input (3d velocity needed by thermal_core) + diagnostic_results->FindResult(&u_g,"u_g"); + + SplitSolutionVectorx(u_g,numberofnodes,3,&vx,&vy,&vz); + femmodel->UpdateInputsFromVector(vx,VxEnum,VertexEnum); + femmodel->UpdateInputsFromVector(vy,VxEnum,VertexEnum); + femmodel->UpdateInputsFromVector(vz,VxEnum,VertexEnum); + + delete diagnostic_results; + } + + /*Compute misfit for this velocity field.*/ + femmodel->UpdateInputsFromConstant(fit[n],FitEnum); + CostFunctionx( &J, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters,analysis_type,sub_analysis_type); + + /*Free ressources:*/ + xfree((void**)&fit); + xfree((void**)&optscal); + xfree((void**)¶m_g_copy); + xfree((void**)&dofset); + xfree((void**)&vx); + xfree((void**)&vy); + xfree((void**)&vz); + + return J; +} Index: /issm/trunk/src/c/solutions/parallel.h =================================================================== --- /issm/trunk/src/c/solutions/parallel.h (revision 3894) +++ /issm/trunk/src/c/solutions/parallel.h (revision 3894) @@ -0,0 +1,63 @@ +/* + * parallel.h: + */ + +#ifndef PARALLEL_H_ +#define PARALLEL_H_ + +#include "../objects/objects.h" +#include "../io/io.h" + +struct OptArgs; +class FemModel; + +/*cores: */ +DataSet* gradjcompute_core(Model* model); +DataSet* diagnostic_core(Model* model); +DataSet* prognostic_core(Model* model); +DataSet* prognostic2_core(Model* model); +DataSet* balancedthickness_core(Model* model); +DataSet* balancedthickness2_core(Model* model); +DataSet* balancedvelocities_core(Model* model); +DataSet* slopecompute_core(Model* model); +DataSet* control_core(Model* model); +DataSet* steadystate_core(Model* model); +DataSet* transient_core(Model* model); +DataSet* transient_core_2d(Model* model); +DataSet* transient_core_3d(Model* model); +DataSet* thermal_core(Model* model); + +/*computational cores: */ +void thermal_core_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem,int analysis_type,int sub_analysis_type); +void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, DataSet* loads, FemModel* fem,int analysis_type,int sub_analysis_type); +void diagnostic_core_linear(Vec* ppg,FemModel* fem,int analysis_type,int sub_analysis_type); +void slope_core(Vec* pslopex,Vec* pslopey,FemModel* fem,int AnalysisEnum); + + + + + +//int GradJOrth(WorkspaceParams* workspaceparams); +void convergence(int* pconverged, Mat K_ff,Vec p_f,Vec u_f,Vec u_f_old,Parameters* parameters); + +int GoldenSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*),FemModel* femmodel); + +int BrentSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*),FemModel* femmodel); + +double objectivefunctionC(double search_scalar,OptArgs* optargs); + +int GradJSearch(double* search_vector,FemModel* femmodel,int step); +//int GradJCheck(WorkspaceParams* workspaceparams,int step,int status); + +//int ParameterUpdate(double* search_vector,int step, WorkspaceParams* workspaceparams,BatchParams* batchparams); +void OutputResults(DataSet* results,char* filename); +void WriteLockFile(char* filename); + +void ControlInitialization(Model* model); +void ControlRestart(Model* model,double* param_g); + +void CreateFemModel(FemModel* femmodel,ConstDataHandle MODEL,int analysis_type,int sub_analysis_type); +//int BatchDebug(Mat* Kgg,Vec* pg,FemModel* femmodel,char* filename); +void ProcessResults(DataSet** pnewresults, DataSet* results,Model* model,int analysis_type); + +#endif Index: /issm/trunk/src/c/solutions/prognostic.cpp =================================================================== --- /issm/trunk/src/c/solutions/prognostic.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/prognostic.cpp (revision 3894) @@ -0,0 +1,133 @@ +/*!\file: prognostic.cpp + * \brief: prognostic solution + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../DataSet/DataSet.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../include/include.h" +#include "../modules.h" +#include "./parallel.h" + +int main(int argc,char* *argv){ + + /*I/O: */ + FILE* fid=NULL; + char* inputfilename=NULL; + char* outputfilename=NULL; + char* lockname=NULL; + int numberofnodes; + bool waitonlock=false; + + Model* model=NULL; + + bool qmu_analysis; + + /*Results: */ + DataSet* results=NULL; + DataSet* processedresults=NULL; + Result* result=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"); + + _printf_("read and create finite element model:\n"); + model->AddFormulation(fid,PrognosticAnalysisEnum); + + /*recover parameters: */ + model->FindParam(&waitonlock,WaitOnLockEnum); + model->FindParam(&qmu_analysis,QmuAnalysisEnum); + + 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 prognostic analysis: */ + _printf_("call computational core:\n"); + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + results=prognostic_core(model); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + + } + else{ + + /*run qmu analysis: */ + _printf_("calling qmu analysis on prognostic core:\n"); + + #ifdef _HAVE_DAKOTA_ + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + Qmux(model,PrognosticAnalysisEnum,NoneAnalysisEnum); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + #else + ISSMERROR(" Dakota not present, cannot do qmu!"); + #endif + } + + + _printf_("process results:\n"); + ProcessResults(&processedresults,results,model,PrognosticAnalysisEnum); + + _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; + + /*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/solutions/prognostic2.cpp =================================================================== --- /issm/trunk/src/c/solutions/prognostic2.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/prognostic2.cpp (revision 3894) @@ -0,0 +1,133 @@ +/*!\file: prognostic2.cpp + * \brief: prognostic2 solution + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../DataSet/DataSet.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../include/include.h" +#include "../modules.h" +#include "./parallel.h" + +int main(int argc,char* *argv){ + + /*I/O: */ + FILE* fid=NULL; + char* inputfilename=NULL; + char* outputfilename=NULL; + char* lockname=NULL; + int numberofnodes; + bool waitonlock=false; + + Model* model=NULL; + + bool qmu_analysis; + + /*Results: */ + DataSet* results=NULL; + DataSet* processedresults=NULL; + Result* result=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"); + + _printf_("read and create finite element model:\n"); + model->AddFormulation(fid,Prognostic2AnalysisEnum); + + /*recover parameters: */ + model->FindParam(&waitonlock,WaitOnLockEnum); + model->FindParam(&qmu_analysis,QmuAnalysisEnum); + model->FindParam(&numberofnodes,NumberOfNodesEnum); + + 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( ); + results=prognostic2_core(model); + 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,Prognostic2AnalysisEnum,NoneAnalysisEnum); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + #else + ISSMERROR(" Dakota not present, cannot do qmu!"); + #endif + } + + _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; + + /*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/solutions/prognostic2_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/prognostic2_core.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/prognostic2_core.cpp (revision 3894) @@ -0,0 +1,77 @@ +/*!\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 "../modules.h" + +DataSet* prognostic2_core(Model* model){ + + extern int my_rank; + + /*output: */ + Result* result=NULL; + DataSet* results=NULL; + + /*intermediary: */ + Vec vx_g=NULL; + Vec vy_g=NULL; + + /*solutions: */ + Vec h_g=NULL; + + /*flags: */ + int verbose=0; + int numberofdofspernode; + int numberofnodes; + int numberofvertices; + int dofs[1]={1}; + + /*fem prognostic model: */ + FemModel* fem_p=NULL; + + //initialize results: + results=new DataSet(ResultsEnum); + + /*recover fem model: */ + fem_p=model->GetFormulation(Prognostic2AnalysisEnum); + + //first recover parameters common to all solutions + model->FindParam(&verbose,VerboseEnum); + model->FindParam(&numberofnodes,NumberOfNodesEnum); + model->FindParam(&numberofvertices,NumberOfVerticesEnum); + model->FindParam(&numberofdofspernode,NumberOfNodesEnum); + + _printf_("depth averaging velocity...\n"); + /*Where is it done?*/ + + _printf_("call computational core:\n"); + diagnostic_core_linear(&h_g,fem_p,Prognostic2AnalysisEnum,NoneAnalysisEnum); + + _printf_("Averaging over vertices:\n"); + FieldAverageOntoVerticesx(&h_g,fem_p->elements,fem_p->nodes,fem_p->vertices,fem_p->loads,fem_p->materials,fem_p->parameters); + + //_printf_("extrude computed thickness on all layers:\n"); + //FieldExtrudex(h_g, fem_p->elements,fem_p->nodes, fem_p->vertices,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); + + /*Add analysis_type to results: */ + result=new Result(results->Size()+1,0,1,"analysis_type",EnumAsString(Prognostic2AnalysisEnum)); + results->AddObject(result); + + /*return: */ + return results; + +} Index: /issm/trunk/src/c/solutions/prognostic_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/prognostic_core.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/prognostic_core.cpp (revision 3894) @@ -0,0 +1,65 @@ +/*!\file: prognostic_core.cpp + * \brief: core of the prognostic solution + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./parallel.h" +#include "../modules.h" + +DataSet* prognostic_core(Model* model){ + + extern int my_rank; + + /*output: */ + DataSet* results=NULL; + Result* result=NULL; + + /*solutions: */ + Vec h_g=NULL; + + /*flags: */ + int verbose=0; + + /*fem prognostic model: */ + FemModel* fem_p=NULL; + + //initialize results + results=new DataSet(ResultsEnum); + + /*recover fem model: */ + fem_p=model->GetFormulation(PrognosticAnalysisEnum); ISSMASSERT(fem_p); + + //first recover parameters common to all solutions + model->FindParam(&verbose,VerboseEnum); + + _printf_("depth averaging velocity...\n"); + DepthAverageInputx(fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,VxEnum); + DepthAverageInputx(fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,VyEnum); + + _printf_("call computational core:\n"); + diagnostic_core_linear(&h_g,fem_p,PrognosticAnalysisEnum,NoneAnalysisEnum); + + _printf_("update inputs:\n"); + fem_p->UpdateInputsFromSolution(h_g,PrognosticAnalysisEnum,NoneAnalysisEnum); + + _printf_("extrude computed thickness on all layers:\n"); + ExtrudeInputx(fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,ThicknessEnum); + + _printf_("extract result from extruded inputs: \n"); + InputToResultx(&result,fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,ThicknessEnum,results->Size()+1,0,1); + results->AddObject(result); + + /*Add analysis_type to results: */ + result=new Result(results->Size()+1,0,1,"analysis_type",EnumAsString(PrognosticAnalysisEnum)); + results->AddObject(result); + + /*Free ressources:*/ + VecFree(&h_g); + + //return: + return results; + +} Index: /issm/trunk/src/c/solutions/slope_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/slope_core.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/slope_core.cpp (revision 3894) @@ -0,0 +1,74 @@ +/*!\file: slope_core.cpp + * \brief: core of the slope solution + */ + +#include "./parallel.h" +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../modules.h" + +void slope_core(Vec* pslopex,Vec* pslopey,FemModel* fem, int AnalysisEnum){ + + /*parameters: */ + int verbose; + int dim; + bool isstokes; + bool ishutter; + int xanalysis; + int yanalysis; + + /*output: */ + Vec slopex=NULL; + Vec slopey=NULL; + + /*Recover some parameters: */ + fem->parameters->FindParam(&verbose,VerboseEnum); + fem->parameters->FindParam(&dim,DimEnum); + fem->parameters->FindParam(&isstokes,IsStokesEnum); + fem->parameters->FindParam(&ishutter,IsHutterEnum); + + if(verbose)_printf_("%s\n","computing slope (x and y derivatives)..."); + + /*Specify type of computations: */ + if(AnalysisEnum==SurfaceAnalysisEnum){ + xanalysis=SurfaceXAnalysisEnum; + yanalysis=SurfaceYAnalysisEnum; + } + else if(AnalysisEnum==BedAnalysisEnum){ + xanalysis=BedXAnalysisEnum; + yanalysis=BedYAnalysisEnum; + } + else ISSMERROR("%s%s%s"," analysis ",EnumAsString(AnalysisEnum)," not supported yet!"); + + + /*Early return possible? */ + if(!ishutter && AnalysisEnum==SurfaceAnalysisEnum){ + /*no need to compute Surface Slope except for Hutter: */ + *pslopex=NULL; + *pslopey=NULL; + return; + } + if(!isstokes && AnalysisEnum==BedAnalysisEnum){ + /*no need to compute Bed Slope except for full Stokes: */ + *pslopex=NULL; + *pslopey=NULL; + return; + } + + + /*Call on core computations: */ + diagnostic_core_linear(&slopex,fem,SlopecomputeAnalysisEnum,xanalysis); + diagnostic_core_linear(&slopey,fem,SlopecomputeAnalysisEnum,yanalysis); + + /*extrude if we are in 3D: */ + if (dim==3){ + if(verbose)_printf_("%s\n","extruding slopes in 3d..."); + FieldExtrudex( slopex, fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,"slopex",0); + FieldExtrudex( slopey, fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,"slopey",0); + } + + /*Assign output pointers:*/ + *pslopex=slopex; + *pslopey=slopey; +} Index: /issm/trunk/src/c/solutions/slopecompute.cpp =================================================================== --- /issm/trunk/src/c/solutions/slopecompute.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/slopecompute.cpp (revision 3894) @@ -0,0 +1,114 @@ +/*!\file: slopecompute.cpp + * \brief: slopecompute solution + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../DataSet/DataSet.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../include/include.h" +#include "./parallel.h" + + +int main(int argc,char* *argv){ + + /*I/O: */ + FILE* fid=NULL; + char* inputfilename=NULL; + char* outputfilename=NULL; + char* lockname=NULL; + int numberofnodes; + bool waitonlock=false; + + Model* model=NULL; + + /*Results: */ + DataSet* results=NULL; + DataSet* processedresults=NULL; + Result* result=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,SlopecomputeAnalysisEnum); + + /*recover parameters: */ + model->FindParam(&waitonlock,WaitOnLockEnum); + + MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime(); + + /*run slopecompute analysis: */ + _printf_("call computational core:\n"); + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + results=slopecompute_core(model); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + + _printf_("process results:\n"); + ProcessResults(&processedresults,results,model,SlopecomputeAnalysisEnum); + + _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; + + /*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/solutions/slopecompute_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/slopecompute_core.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/slopecompute_core.cpp (revision 3894) @@ -0,0 +1,67 @@ +/*!\file: slopecompute_core.cpp + * \brief: core of the slopecompute solution + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./parallel.h" +#include "../modules.h" + +DataSet* slopecompute_core(Model* model){ + + extern int my_rank; + + /*output: */ + DataSet* results=NULL; + Result* result=NULL; + + /*solutions: */ + Vec sx_g=NULL; + Vec sy_g=NULL; + + /*flags: */ + int verbose=0; + int numberofdofspernode; + int numberofnodes; + int dofs[2]={1,1}; + + /*fem slopecompute model: */ + FemModel* fem_sl=NULL; + + //initialize results + results=new DataSet(ResultsEnum); + + /*recover fem model: */ + fem_sl=model->GetFormulation(SlopecomputeAnalysisEnum); + + //first recover parameters common to all solutions + model->FindParam(&verbose,VerboseEnum); + model->FindParam(&numberofnodes,NumberOfNodesEnum); + model->FindParam(&numberofdofspernode,NumberOfDofsPerNodeEnum); + + _printf_("call computational core:\n"); + diagnostic_core_linear(&sx_g,fem_sl,SlopecomputeAnalysisEnum,SurfaceXAnalysisEnum); + diagnostic_core_linear(&sy_g,fem_sl,SlopecomputeAnalysisEnum,SurfaceYAnalysisEnum); + + _printf_("extrude computed slope on all layers:\n"); + FieldExtrudex( sx_g, fem_sl->elements,fem_sl->nodes, fem_sl->vertices,fem_sl->loads, fem_sl->materials,fem_sl->parameters,"slopex",0); + FieldExtrudex( sy_g, fem_sl->elements,fem_sl->nodes, fem_sl->vertices,fem_sl->loads, fem_sl->materials,fem_sl->parameters,"slopey",0); + + /*Plug results into output dataset: */ + result=new Result(results->Size()+1,0,1,"sx_g",sx_g); + results->AddObject(result); + result=new Result(results->Size()+1,0,1,"sy_g",sy_g); + results->AddObject(result); + result=new Result(results->Size()+1,0,1,"analysis_type",EnumAsString(SlopecomputeAnalysisEnum)); + results->AddObject(result); + + + /*Free ressources:*/ + VecFree(&sx_g); + VecFree(&sy_g); + + //return: + return results; +} Index: /issm/trunk/src/c/solutions/steadystate.cpp =================================================================== --- /issm/trunk/src/c/solutions/steadystate.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/steadystate.cpp (revision 3894) @@ -0,0 +1,182 @@ +/*!\file: steadystate.cpp + * \brief: steadystate solution + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../DataSet/DataSet.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../include/include.h" +#include "../modules.h" +#include "./parallel.h" + +int main(int argc,char* *argv){ + + /*I/O: */ + FILE* fid=NULL; + char* inputfilename=NULL; + char* outputfilename=NULL; + char* lockname=NULL; + int numberofnodes; + bool qmu_analysis=false; + bool control_analysis=false; + char* control_type=NULL; + + /*Model: */ + Model* model=NULL; + FemModel* fem_dh=NULL; + FemModel* fem_ds=NULL; + + /*Results: */ + DataSet* results=NULL; + DataSet* processed_results=NULL; + Result* result=NULL; + + bool waitonlock=false; + + double* u_g_initial=NULL; + double* p_g_initial=NULL; + double* u_g_obs=NULL; + double* weights=NULL; + double dt; + BoolParam* 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]; + + /*Open handle to data on disk: */ + fid=pfopen(inputfilename,"rb"); + + /*Initialize model structure: */ + MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime(); + model=new Model(); + + _printf_("read and create finite element model:\n"); + + _printf_("\n reading diagnostic horiz model data:\n"); + model->AddFormulation(fid,DiagnosticAnalysisEnum,HorizAnalysisEnum); + + _printf_("\n reading diagnostic vert model data:\n"); + model->AddFormulation(fid,DiagnosticAnalysisEnum,VertAnalysisEnum); + + _printf_("\n reading diagnostic stokes model data:\n"); + model->AddFormulation(fid,DiagnosticAnalysisEnum,StokesAnalysisEnum); + + _printf_("\n reading diagnostic hutter model data:\n"); + model->AddFormulation(fid,DiagnosticAnalysisEnum,HutterAnalysisEnum); + + _printf_("\n reading surface and bed slope computation model data:\n"); + model->AddFormulation(fid,SlopecomputeAnalysisEnum); + + _printf_("\n read and create thermal finite element model:\n"); + model->AddFormulation(fid,ThermalAnalysisEnum); + _printf_("\n read and create melting finite element model:\n"); + model->AddFormulation(fid,MeltingAnalysisEnum); + + /*recover parameters: */ + model->FindParam(&waitonlock,WaitOnLockEnum); + model->FindParam(&control_analysis,ControlAnalysisEnum); + model->FindParam(&qmu_analysis,QmuAnalysisEnum); + + MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime(); + + /*are we running the solution sequence, or a qmu wrapper around it? : */ + if(!qmu_analysis){ + + if(!control_analysis){ + /*run diagnostic analysis: */ + _printf_("call computational core:\n"); + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + results=steadystate_core(model); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + + _printf_("process results:\n"); + ProcessResults(&processed_results,results,model,SteadystateAnalysisEnum); + } + else{ + /*change control_steady to 1 to know we are doing steadystate*/ + fem_dh=model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum); + fem_ds=model->GetFormulation(DiagnosticAnalysisEnum,StokesAnalysisEnum); + + param=(BoolParam*)fem_dh->parameters->FindParamObject(ControlSteadyEnum); param->value=true; + param=(BoolParam*)fem_ds->parameters->FindParamObject(ControlSteadyEnum); param->value=true; + + /*run control analysis: */ + _printf_("call computational core:\n"); + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + results=control_core(model); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + + _printf_("process results:\n"); + ProcessResults(&processed_results,results,model,ControlAnalysisEnum); + } + + _printf_("write results to disk:\n"); + OutputResults(processed_results,outputfilename); + } + else{ + /*run qmu analysis: */ + _printf_("calling qmu analysis on steadystate core:\n"); + + #ifdef _HAVE_DAKOTA_ + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + Qmux(model,SteadystateAnalysisEnum,NoneAnalysisEnum); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + #else + ISSMERROR(" Dakota not present, cannot do qmu!"); + #endif + } + + _printf_("write lock file:\n"); + if (waitonlock){ + WriteLockFile(lockname); + } + + /*Free ressources */ + xfree((void**)&u_g_initial); + xfree((void**)&u_g_obs); + xfree((void**)&p_g_initial); + delete model; + delete results; + delete processed_results; + + /*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/solutions/steadystate_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/steadystate_core.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/steadystate_core.cpp (revision 3894) @@ -0,0 +1,149 @@ +/*!\file: steadystate_core.cpp + * \brief: core of the steadystate solution + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./parallel.h" +#include "../modules.h" + +DataSet* steadystate_core(Model* model){ + + extern int my_rank; + + /*fem models: */ + FemModel* fem_dh=NULL; + FemModel* fem_dv=NULL; + FemModel* fem_dhu=NULL; + FemModel* fem_ds=NULL; + FemModel* fem_sl=NULL; + FemModel* fem_t=NULL; + FemModel* fem_m=NULL; + + /*output: */ + Result* result=NULL; + DataSet* results=NULL; + DataSet* results_thermal=NULL; + DataSet* results_diagnostic=NULL; + + /*solutions: */ + Vec u_g=NULL; + Vec old_u_g=NULL; + Vec t_g=NULL; + Vec t_g_average=NULL; + Vec old_t_g=NULL; + Vec p_g=NULL; + Vec m_g=NULL; + Vec du_g=NULL; + Vec dt_g=NULL; + double ndu,nu; + double normdt,normt; + double eps_rel; + + /*flags: */ + int verbose=0; + int isstokes=0; + int numberofnodes; + int ndof; + int converged; + int step; + + /*recover fem models: */ + fem_dh=model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum); + fem_dv=model->GetFormulation(DiagnosticAnalysisEnum,VertAnalysisEnum); + fem_ds=model->GetFormulation(DiagnosticAnalysisEnum,StokesAnalysisEnum); + fem_dhu=model->GetFormulation(DiagnosticAnalysisEnum,HutterAnalysisEnum); + fem_sl=model->GetFormulation(SlopecomputeAnalysisEnum); + fem_t=model->GetFormulation(ThermalAnalysisEnum); + fem_m=model->GetFormulation(MeltingAnalysisEnum); + + + //first recover parameters common to all solutions + model->FindParam(&verbose,VerboseEnum); + model->FindParam(&numberofnodes,NumberOfNodesEnum); + model->FindParam(&eps_rel,EpsResEnum); + model->FindParam(&isstokes,IsStokesEnum); + + //initialize: + converged=0; + step=1; + + if (isstokes)ndof=4; + else ndof=3; + + for(;;){ + + if(verbose)_printf_("%s%i\n"," computing temperature and velocity for step: ",step); + + //first compute temperature at steady state. + results_thermal=thermal_core(model); + + //get t_g and m_g; + VecFree(&t_g);results_thermal->FindResult(&t_g,"t_g"); + VecFree(&m_g);results_thermal->FindResult(&m_g,"m_g"); + delete results_thermal; + + //Add temperature to inputs. + //compute depth averaged temperature and add to inputs + VecDuplicatePatch(&t_g_average,t_g); + FieldDepthAveragex( t_g_average, fem_t->elements,fem_t->nodes, fem_t->vertices,fem_t->loads, fem_t->materials,fem_t->parameters,"temperature"); + model->UpdateInputsFromVector(t_g_average,TemperatureAverageEnum,VertexEnum); + model->UpdateInputsFromVector(t_g,TemperatureEnum,VertexEnum); + VecFree(&t_g_average); //not needed anymore + + //now compute diagnostic velocity using the steady state temperature. + results_diagnostic=diagnostic_core(model); + + //get p_g and u_g + VecFree(&u_g);results_diagnostic->FindResult(&u_g,"u_g"); + VecFree(&p_g);results_diagnostic->FindResult(&p_g,"p_g"); + delete results_diagnostic; + + //convergence? + if(step>1){ + VecDuplicatePatch(&du_g,old_u_g);VecAYPX(du_g,-1.0,u_g); + VecNorm(du_g,NORM_2,&ndu); VecNorm(old_u_g,NORM_2,&nu); VecFree(&du_g); + + VecDuplicatePatch(&dt_g,old_t_g); VecAYPX(dt_g,-1.0,t_g); + VecNorm(dt_g,NORM_2,&normdt); VecNorm(old_t_g,NORM_2,&normt);VecFree(&dt_g); + + if (verbose) _printf_("%-60s%g\n %s%g\n %s%g%s\n", + " relative convergence criterion: velocity -> norm(du)/norm(u)= ",ndu/nu*100," temperature -> norm(dt)/norm(t)=",normdt/normt*100," eps_rel: ",eps_rel*100," %"); + + if ((ndu/nu<=eps_rel) && (normdt/normt<=eps_rel)) converged=1; + else converged=0; + } + else{ + converged=0; + } + + VecFree(&old_u_g);VecDuplicatePatch(&old_u_g,u_g); + VecFree(&old_t_g);VecDuplicatePatch(&old_t_g,t_g); + + step++; + if (converged)break; + } + + /*Plug results into output dataset: */ + result=new Result(results->Size()+1,0,1,"u_g",u_g); + results->AddObject(result); + result=new Result(results->Size()+1,0,1,"p_g",p_g); + results->AddObject(result); + result=new Result(results->Size()+1,0,1,"t_g",t_g); + results->AddObject(result); + result=new Result(results->Size()+1,0,1,"m_g",m_g); + results->AddObject(result); + result=new Result(results->Size()+1,0,1,"analysis_type",EnumAsString(SteadyAnalysisEnum)); + results->AddObject(result); + + + /*Free ressource*/ + VecFree(&old_u_g); + VecFree(&old_t_g); + VecFree(&u_g); + VecFree(&p_g); + VecFree(&t_g); + VecFree(&m_g); +} Index: /issm/trunk/src/c/solutions/thermal.cpp =================================================================== --- /issm/trunk/src/c/solutions/thermal.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/thermal.cpp (revision 3894) @@ -0,0 +1,135 @@ +/*!\file: thermal.cpp + * \brief: thermal solution + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../DataSet/DataSet.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../include/include.h" +#include "../modules.h" +#include "./parallel.h" + +int main(int argc,char* *argv){ + + int i,n; + + /*I/O: */ + FILE* fid=NULL; + char* inputfilename=NULL; + char* outputfilename=NULL; + char* lockname=NULL; + bool qmu_analysis=false; + int numberofnodes; + bool waitonlock=false; + + /*Model: */ + Model* model=NULL; + + /*Results: */ + DataSet* results=NULL; + DataSet* processed_results=NULL; + Result* result=NULL; + + Param* param=NULL; + double dt; + double yts; + + /*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); + + inputfilename=argv[2]; + outputfilename=argv[3]; + lockname=argv[4]; + + /*Open handle to data on disk: */ + fid=pfopen(inputfilename,"rb"); + + /*Initialize model structure: */ + MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime(); + model=new Model(); + + _printf_("read and create thermal finite element model:\n"); + model->AddFormulation(fid,ThermalAnalysisEnum); + _printf_("read and create melting finite element model:\n"); + model->AddFormulation(fid,MeltingAnalysisEnum); + + /*recover parameters: */ + model->FindParam(&waitonlock,WaitOnLockEnum); + model->FindParam(&qmu_analysis,QmuAnalysisEnum); + + 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 thermal analysis: */ + _printf_("call computational core:\n"); + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + results=thermal_core(model); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + + _printf_("process results:\n"); + ProcessResults(&processed_results,results,model,ThermalAnalysisEnum); + + _printf_("write results to disk:\n"); + OutputResults(processed_results,outputfilename); + } + else{ + /*run qmu analysis: */ + _printf_("calling qmu analysis on thermal core:\n"); + + #ifdef _HAVE_DAKOTA_ + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + Qmux(model,ThermalAnalysisEnum,NoneAnalysisEnum); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + #else + ISSMERROR(" Dakota not present, cannot do qmu!"); + #endif + } + + if (waitonlock>0){ + _printf_("write lock file:\n"); + WriteLockFile(lockname); + } + + /*Free ressources */ + delete model; + delete results; + delete processed_results; + + /*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/solutions/thermal_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/thermal_core.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/thermal_core.cpp (revision 3894) @@ -0,0 +1,139 @@ +/*!\file: thermal_core.cpp + * \brief: core of the thermal solution + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./parallel.h" +#include "../modules.h" + +DataSet* thermal_core(Model* model){ + + extern int my_rank; + int i; + + /*fem models: */ + FemModel* fem_t=NULL; + FemModel* fem_m=NULL; + + /*output: */ + DataSet* results=NULL; + Result* result=NULL; + + /*solutions vectors: */ + Vec* t_g=NULL; + Vec* m_g=NULL; + double* time=NULL; + double* t_g_serial=NULL; + double* m_g_serial=NULL; + + /*intermediary*/ + Vec Vect=NULL; + Vec Vecm=NULL; + + /*flags: */ + int verbose=0; + int numberofdofspernode; + int numberofnodes; + int nsteps; + double ndt; + double dt; + + int sub_analysis_type; + double melting_offset; + + //initialize results + results=new DataSet(ResultsEnum); + + /*recover fem models: */ + fem_t=model->GetFormulation(ThermalAnalysisEnum); + fem_m=model->GetFormulation(MeltingAnalysisEnum); + + //first recover parameters common to all solutions + fem_t->FindParam(&numberofnodes,NumberOfNodesEnum); + fem_t->FindParam(&sub_analysis_type,SubAnalysisTypeEnum); + fem_t->FindParam(&verbose,VerboseEnum); + fem_t->FindParam(&ndt,NdtEnum); + fem_t->FindParam(&dt,DtEnum); + + if(dt==0){ + + time=(double*)xmalloc(sizeof(double)); + time[0]=0; + + /*allocate t_g and m_g arrays: */ + t_g=(Vec*)xmalloc(sizeof(Vec)); + m_g=(Vec*)xmalloc(sizeof(Vec)); + + if(verbose)_printf_("computing temperatures:\n"); + thermal_core_nonlinear(&t_g[0],&melting_offset,fem_t,ThermalAnalysisEnum,NoneAnalysisEnum); + model->UpdateInputsFromVector(t_g[0],TemperatureEnum,VertexEnum); + model->UpdateInputsFromConstant(melting_offset,MeltingOffsetEnum); + + if(verbose)_printf_("computing melting:\n"); + diagnostic_core_linear(&m_g[0],fem_m,MeltingAnalysisEnum,NoneAnalysisEnum); + } + else{ + + nsteps=(int)(ndt/dt); + time=(double*)xmalloc((nsteps)*sizeof(double)); + + /*allocate t_g and m_g arrays: */ + t_g=(Vec*)xmalloc((nsteps)*sizeof(Vec)); + m_g=(Vec*)xmalloc((nsteps)*sizeof(Vec)); + + for(i=0;iUpdateInputsFromVector(t_g[i],TemperatureEnum,VertexEnum); + model->UpdateInputsFromConstant(melting_offset,MeltingOffsetEnum); + diagnostic_core_linear(&m_g[i],fem_m,MeltingAnalysisEnum,NoneAnalysisEnum); + } + } + + /*Plug results into output dataset: */ + if(dt==0){ + result=new Result(results->Size()+1,0,1,"t_g",t_g[0]); + results->AddObject(result); + + result=new Result(results->Size()+1,0,1,"m_g",m_g[0]); + results->AddObject(result); + + /*free ressource*/ + Vect=t_g[0]; Vecm=m_g[0]; + VecFree(&Vect); VecFree(&Vecm); + } + else{ + for(i=0;iSize()+1,time[i],i+1,"t_g",t_g[i]); + results->AddObject(result); + + result=new Result(results->Size()+1,time[i],i+1,"m_g",m_g[i]); + results->AddObject(result); + + /*free ressource*/ + Vect=t_g[i]; Vecm=m_g[i]; + VecFree(&Vect); VecFree(&Vecm); + } + } + /*Add analysis_type to results: */ + result=new Result(results->Size()+1,0,1,"analysis_type",EnumAsString(ThermalAnalysisEnum)); + results->AddObject(result); + + /*free ressource*/ + xfree((void**)&t_g); + xfree((void**)&m_g); + xfree((void**)&time); + + //return: + return results; + + +} Index: /issm/trunk/src/c/solutions/thermal_core_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/solutions/thermal_core_nonlinear.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/thermal_core_nonlinear.cpp (revision 3894) @@ -0,0 +1,133 @@ +/*!\file: thermal_core_nonlinear.cpp + * \brief: core of the thermal solution + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../modules.h" + +void thermal_core_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem,int analysis_type,int sub_analysis_type){ + + /*solution : */ + Vec tg=NULL; + Vec tf=NULL; + Vec tf_old=NULL; + double melting_offset; + + /*intermediary: */ + Mat Kgg=NULL; + Mat Kgg_nopenalty=NULL; + Mat Kff=NULL; + Mat Kfs=NULL; + Vec pg=NULL; + Vec pg_nopenalty=NULL; + Vec pf=NULL; + + int converged; + int constraints_converged; + int num_unstable_constraints; + int count; + int numberofnodes; + int min_thermal_constraints; + bool reset_penalties; + + /*parameters:*/ + int kflag,pflag; + char* solver_string=NULL; + int verbose=0; + bool lowmem=0; + + /*Recover parameters: */ + kflag=1; pflag=1; + + fem->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + fem->parameters->FindParam(&solver_string,SolverStringEnum); + fem->parameters->FindParam(&verbose,VerboseEnum); + fem->parameters->FindParam(&lowmem,LowmemEnum); + fem->parameters->FindParam(&min_thermal_constraints,MinThermalConstraintsEnum); + + count=1; + converged=0; + + for(;;){ + + if(verbose)_printf_("%s\n","starting direct shooting method"); + + if(count==1) reset_penalties=1; else reset_penalties=0; + fem->UpdateInputsFromConstant(reset_penalties,ResetPenaltiesEnum); + + //*Generate system matrices + if (!lowmem){ + + /*Compute Kgg_nopenalty and pg_nopenalty once for all: */ + if (count==1){ + SystemMatricesx(&Kgg_nopenalty, &pg_nopenalty,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); + } + + /*Copy K_gg_nopenalty into Kgg, same for pg: */ + MatDuplicate(Kgg_nopenalty,MAT_COPY_VALUES,&Kgg); + VecDuplicatePatch(&pg,pg_nopenalty); + + //apply penalties each time + PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); + } + else{ + SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); + //apply penalties + PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); + } + + /*!Reduce matrix from g to f size:*/ + Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets); + + /*Free ressources: */ + MatFree(&Kgg); + + if (verbose) _printf_(" reducing load from g to f set\n"); + /*!Reduce load from g to f size: */ + Reduceloadfromgtofx(&pf, pg, fem->Gmn, Kfs, fem->ys, fem->nodesets); + + //no need for pg and Kfs anymore + VecFree(&pg); + MatFree(&Kfs); + + /*Solve: */ + if(verbose)_printf_("%s\n","solving"); + VecFree(&tf); + Solverx(&tf, Kff, pf,tf_old, solver_string); + VecFree(&tf_old); VecDuplicatePatch(&tf_old,tf); + + //no need for Kff and pf anymore + MatFree(&Kff);VecFree(&pf);VecFree(&tg); + + if (verbose) _printf_(" merging solution from f to g set\n"); + //Merge back to g set + Mergesolutionfromftogx(&tg, tf,fem->Gmn,fem->ys,fem->nodesets); + + //Deal with penalty loads + if (verbose) _printf_(" penalty constraints\n"); + PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,analysis_type,sub_analysis_type); + + fem->UpdateInputsFromVector(tg,TemperatureEnum,VertexEnum); + + if (!converged){ + if(verbose)_printf_("%s%i\n"," #unstable constraints = ",num_unstable_constraints); + if (num_unstable_constraints <= min_thermal_constraints)converged=1; + } + count++; + + if(converged==1)break; + } + + /*Free ressources: */ + MatFree(&Kgg_nopenalty); + VecFree(&pg_nopenalty); + VecFree(&tf); + VecFree(&tf_old); + delete solver_string; + + /*Assign output pointers: */ + *ptg=tg; + *pmelting_offset=melting_offset; +} Index: /issm/trunk/src/c/solutions/transient.cpp =================================================================== --- /issm/trunk/src/c/solutions/transient.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/transient.cpp (revision 3894) @@ -0,0 +1,156 @@ +/*!\file: transient.cpp + * \brief: transient solution + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../DataSet/DataSet.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../include/include.h" +#include "../modules.h" +#include "./parallel.h" + +int main(int argc,char* *argv){ + + /*I/O: */ + FILE* fid=NULL; + char* inputfilename=NULL; + char* outputfilename=NULL; + char* lockname=NULL; + int numberofnodes; + bool qmu_analysis=false; + bool waitonlock=false; + + /*Model: */ + Model* model=NULL; + int dim=-1; + + /*Results: */ + DataSet* results=NULL; + DataSet* processed_results=NULL; + Result* result=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]; + + /*Open handle to data on disk: */ + fid=pfopen(inputfilename,"rb"); + + /*Initialize model structure: */ + MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime(); + model=new Model(); + + _printf_("read and create finite element model:\n"); + _printf_("\n reading diagnostic horiz model data:\n"); + model->AddFormulation(fid,DiagnosticAnalysisEnum,HorizAnalysisEnum); + + _printf_("\n reading diagnostic vert model data:\n"); + model->AddFormulation(fid,DiagnosticAnalysisEnum,VertAnalysisEnum); + + _printf_("\n reading diagnostic stokes model data:\n"); + model->AddFormulation(fid,DiagnosticAnalysisEnum,StokesAnalysisEnum); + + _printf_("\n reading diagnostic hutter model data:\n"); + model->AddFormulation(fid,DiagnosticAnalysisEnum,HutterAnalysisEnum); + + _printf_("\n reading surface and bed slope computation model data:\n"); + model->AddFormulation(fid,SlopecomputeAnalysisEnum); + + _printf_("\n reading prognositc model data:\n"); + model->AddFormulation(fid,PrognosticAnalysisEnum); + + /*Do we run in 3d?, in which case we need thermal and melting also:*/ + model->FindParam(&dim,DimEnum); + if(dim==3){ + _printf_("read and create thermal finite element model:\n"); + model->AddFormulation(fid,ThermalAnalysisEnum,TransientAnalysisEnum); + _printf_("read and create melting finite element model:\n"); + model->AddFormulation(fid,MeltingAnalysisEnum,TransientAnalysisEnum); + } + + /*recover parameters: */ + model->FindParam(&waitonlock,WaitOnLockEnum); + model->FindParam(&qmu_analysis,QmuAnalysisEnum); + + MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime(); + + /*are we running the solution sequence, or a qmu wrapper around it? : */ + if(!qmu_analysis){ + + /*run diagnostic analysis: */ + _printf_("call computational core:\n"); + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + results=transient_core(model); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + + _printf_("process results:\n"); + ProcessResults(&processed_results,results,model,TransientAnalysisEnum); + + _printf_("write results to disk:\n"); + OutputResults(processed_results,outputfilename); + } + else{ + /*run qmu analysis: */ + _printf_("calling qmu analysis on transient core:\n"); + + #ifdef _HAVE_DAKOTA_ + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + Qmux(model,TransientAnalysisEnum,NoneAnalysisEnum); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + #else + ISSMERROR(" Dakota not present, cannot do qmu!"); + #endif + } + + if (waitonlock>0){ + _printf_("write lock file:\n"); + WriteLockFile(lockname); + } + + /*Free ressources:*/ + delete results; + delete processed_results; + delete model; + + /*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/solutions/transient_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/transient_core.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/transient_core.cpp (revision 3894) @@ -0,0 +1,29 @@ +/*!\file: transient_core.cpp + * \brief: core of the transient solution + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./parallel.h" +#include "../modules.h" + +DataSet* transient_core(Model* model){ + + int dim=-1; + DataSet* results=NULL; + + //first recover parameters common to all solutions + model->FindParam(&dim,DimEnum); + + //branch out + if(dim==2){ + results=transient_core_2d(model); + } + else if(dim==3){ + results=transient_core_3d(model); + } + else ISSMERROR("%s%i%s"," dimension ",dim," not supported yet!"); + +} Index: /issm/trunk/src/c/solutions/transient_core_2d.cpp =================================================================== --- /issm/trunk/src/c/solutions/transient_core_2d.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/transient_core_2d.cpp (revision 3894) @@ -0,0 +1,156 @@ +/*!\file: transient_2d_core.cpp + * \brief: core of the transient_2d solution + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./parallel.h" +#include "../modules.h" + +DataSet* transient_core_2d(Model* model){ + + extern int my_rank; + + /*fem models: */ + FemModel* fem_dh=NULL; + FemModel* fem_dv=NULL; + FemModel* fem_dhu=NULL; + FemModel* fem_ds=NULL; + FemModel* fem_sl=NULL; + FemModel* fem_p=NULL; + + /*output: */ + DataSet* results=NULL; + Result* result=NULL; + + /*solutions: */ + int step,i; + double time; + + Vec u_g=NULL; + Vec p_g=NULL; + Vec h_g=NULL; + Vec new_h_g=NULL; + Vec h_g_intermediary=NULL; + Vec s_g=NULL; + Vec new_s_g=NULL; + Vec b_g=NULL; + Vec new_b_g=NULL; + double* vx=NULL; + double* vy=NULL; + + /*intermediary: */ + double finaltime; + double dt; + DataSet* diagnostic_results=NULL; + DataSet* prognostic_results=NULL; + + + /*flags: */ + int verbose=0; + int numberofnodes; + + /*dof recovery: */ + int dof01[2]={0,1}; + int dof2[1]={2}; + int dof012[3]={0,1,2}; + int dof3[1]={3}; + + //initialize results + results=new DataSet(ResultsEnum); + + /*recover fem models: */ + fem_dh=model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum); + fem_dv=model->GetFormulation(DiagnosticAnalysisEnum,VertAnalysisEnum); + fem_ds=model->GetFormulation(DiagnosticAnalysisEnum,StokesAnalysisEnum); + fem_dhu=model->GetFormulation(DiagnosticAnalysisEnum,HutterAnalysisEnum); + fem_sl=model->GetFormulation(SlopecomputeAnalysisEnum); + fem_p=model->GetFormulation(PrognosticAnalysisEnum); + + + //first recover parameters common to all solutions + model->FindParam(&verbose,VerboseEnum); + model->FindParam(&finaltime,NdtEnum); + model->FindParam(&dt,DtEnum); + model->FindParam(&numberofnodes,NumberOfNodesEnum); + + /*initialize: */ + step=1; + time=0; + + while(timeFindResult(&u_g,"u_g"); + VecFree(&p_g); diagnostic_results->FindResult(&p_g,"p_g"); + delete diagnostic_results; + + //compute new thickness + _printf_("%s\n","computing new thickness"); + + SplitSolutionVectorx(u_g,numberofnodes,2,&vx,&vy); + model->UpdateInputsFromVector(vx,VxEnum,VertexEnum); + model->UpdateInputsFromVector(vy,VyEnum,VertexEnum); + xfree((void**)&vx); xfree((void**)&vy); + + prognostic_results= prognostic_core(model); + + //extract h_g prognostic_results, and erase prognostic_results; + prognostic_results->FindResult(&h_g_intermediary,"h_g"); + delete prognostic_results; + + //update surface and bed using the new thickness + _printf_(" updating geometry\n"); + UpdateGeometryx(&new_h_g,&new_b_g,&new_s_g, fem_p->elements, fem_p->nodes,fem_p->vertices,fem_p->loads, fem_p->materials, fem_p->parameters,h_g_intermediary,b_g,s_g); + + VecFree(&h_g);h_g=new_h_g; + VecFree(&s_g);s_g=new_s_g; + VecFree(&b_g);b_g=new_b_g; + VecFree(&h_g_intermediary); + + //plug into results. + result=new Result(results->Size()+1,time,step,"u_g",u_g); results->AddObject(result); + result=new Result(results->Size()+1,time,step,"p_g",p_g); results->AddObject(result); + result=new Result(results->Size()+1,time,step,"h_g",h_g); results->AddObject(result); + result=new Result(results->Size()+1,time,step,"s_g",s_g); results->AddObject(result); + result=new Result(results->Size()+1,time,step,"b_g",b_g); results->AddObject(result); + + //update inputs + model->UpdateInputsFromVector(h_g,ThicknessEnum,VertexEnum); + model->UpdateInputsFromVector(s_g,SurfaceEnum,VertexEnum); + model->UpdateInputsFromVector(b_g,BedEnum,VertexEnum); + + SplitSolutionVectorx(u_g,numberofnodes,2,&vx,&vy); + model->UpdateInputsFromVector(vx,VxEnum,VertexEnum); + model->UpdateInputsFromVector(vy,VyEnum,VertexEnum); + + } + + /*Add analysis_type to results: */ + result=new Result(results->Size()+1,0,1,"analysis_type",EnumAsString(TransientAnalysisEnum)); + results->AddObject(result); + + /*Free ressources:*/ + VecFree(&u_g); + VecFree(&p_g); + VecFree(&h_g); + VecFree(&s_g); + VecFree(&b_g); + xfree((void**)&vx); + xfree((void**)&vy); + + //return: + return results; + +} Index: /issm/trunk/src/c/solutions/transient_core_3d.cpp =================================================================== --- /issm/trunk/src/c/solutions/transient_core_3d.cpp (revision 3894) +++ /issm/trunk/src/c/solutions/transient_core_3d.cpp (revision 3894) @@ -0,0 +1,211 @@ +/*!\file: transient_3d_core.cpp + * \brief: core of the transient_3d solution + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./parallel.h" +#include "../modules.h" + +DataSet* transient_core_3d(Model* model){ + + extern int my_rank; + + /*fem models: */ + FemModel* fem_dh=NULL; + FemModel* fem_dv=NULL; + FemModel* fem_dhu=NULL; + FemModel* fem_ds=NULL; + FemModel* fem_sl=NULL; + FemModel* fem_p=NULL; + FemModel* fem_t=NULL; + FemModel* fem_m=NULL; + + /*output: */ + DataSet* results=NULL; + Result* result=NULL; + + /*solutions: */ + int step; + double time; + + double* vx=NULL; + double* vy=NULL; + double* vz=NULL; + + Vec u_g=NULL; + Vec p_g=NULL; + Vec h_g=NULL; + Vec h_g_intermediary=NULL; + Vec s_g=NULL; + Vec b_g=NULL; + Vec t_g=NULL; + Vec t_g_average=NULL; + double melting_offset; + Vec m_g=NULL; + + /*intermediary: */ + double finaltime; + double dt; + DataSet* diagnostic_results=NULL; + DataSet* prognostic_results=NULL; + + + /*flags: */ + int verbose=0; + int numberofnodes; + + /*dof recovery: */ + int dof01[2]={0,1}; + int dof2[1]={2}; + int dof012[3]={0,1,2}; + int dof3[1]={3}; + + //initialize results + results=new DataSet(ResultsEnum); + + /*recover fem models: */ + fem_dh=model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum); + fem_dv=model->GetFormulation(DiagnosticAnalysisEnum,VertAnalysisEnum); + fem_ds=model->GetFormulation(DiagnosticAnalysisEnum,StokesAnalysisEnum); + fem_dhu=model->GetFormulation(DiagnosticAnalysisEnum,HutterAnalysisEnum); + fem_sl=model->GetFormulation(SlopecomputeAnalysisEnum); + fem_p=model->GetFormulation(PrognosticAnalysisEnum); + fem_t=model->GetFormulation(ThermalAnalysisEnum); + fem_m=model->GetFormulation(MeltingAnalysisEnum); + + + //first recover parameters common to all solutions + model->FindParam(&verbose,VerboseEnum); + model->FindParam(&finaltime,NdtEnum); + model->FindParam(&dt,DtEnum); + model->FindParam(&numberofnodes,NumberOfDofsPerNodeEnum); + + /*initialize: */ + step=1; + time=0; + + while(timeUpdateInputsFromVector(t_g,TemperatureEnum,VertexEnum); + model->UpdateInputsFromConstant(melting_offset,MeltingOffsetEnum); + diagnostic_core_linear(&m_g,fem_m,MeltingAnalysisEnum,TransientAnalysisEnum); + + //Compute depth averaged temperature and add to inputs + if(verbose)_printf_("%s\n","computing depth average temperature"); + VecDuplicatePatch(&t_g_average,t_g); + FieldDepthAveragex( t_g_average, fem_t->elements,fem_t->nodes, fem_t->vertices,fem_t->loads, fem_t->materials,fem_t->parameters,"temperature"); + model->UpdateInputsFromVector(t_g_average,TemperatureAverageEnum,VertexEnum); + VecFree(&t_g_average); //not needed anymore + + //Deal with velocities. + diagnostic_results=diagnostic_core(model); + + //extract u_g and p_g from diagnostic_results, and erase diagnostic_results; + diagnostic_results->FindResult(&u_g,"u_g"); + diagnostic_results->FindResult(&p_g,"p_g"); + delete diagnostic_results; + + //compute new thickness + if(verbose)_printf_("%s\n","computing new thickness"); + + SplitSolutionVectorx(u_g,numberofnodes,3,&vx,&vy,&vz); + model->UpdateInputsFromVector(vx,VxEnum,VertexEnum); + model->UpdateInputsFromVector(vy,VyEnum,VertexEnum); + model->UpdateInputsFromVector(vz,VzEnum,VertexEnum); + xfree((void**)&vx); xfree((void**)&vy); xfree((void**)&vz); + + prognostic_results=prognostic_core(model); + + //extract h_g prognostic_results, and erase prognostic_results; + prognostic_results->FindResult(&h_g_intermediary,"h_g"); + delete prognostic_results; + + //update surface and bed using the new thickness + if(verbose)_printf_(" updating geometry\n"); + UpdateGeometryx(&h_g,&b_g,&s_g, fem_p->elements, fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials, fem_p->parameters,h_g_intermediary,b_g,s_g); + VecFree(&h_g_intermediary); + + if(verbose)_printf_("%s\n","updating node positions"); + UpdateVertexPositionsx( fem_dh ->vertices,h_g,b_g); + UpdateVertexPositionsx( fem_dv ->vertices,h_g,b_g); + UpdateVertexPositionsx( fem_dhu->vertices,h_g,b_g); + UpdateVertexPositionsx( fem_ds ->vertices,h_g,b_g); + UpdateVertexPositionsx( fem_sl ->vertices,h_g,b_g); + UpdateVertexPositionsx( fem_p ->vertices,h_g,b_g); + UpdateVertexPositionsx( fem_t ->vertices,h_g,b_g); + UpdateVertexPositionsx( fem_m ->vertices,h_g,b_g); + + //plug into results. + result=new Result(results->Size()+1,time,step,"u_g",u_g); results->AddObject(result); + result=new Result(results->Size()+1,time,step,"p_g",p_g); results->AddObject(result); + result=new Result(results->Size()+1,time,step,"h_g",h_g); results->AddObject(result); + result=new Result(results->Size()+1,time,step,"s_g",s_g); results->AddObject(result); + result=new Result(results->Size()+1,time,step,"b_g",b_g); results->AddObject(result); + result=new Result(results->Size()+1,time,step,"t_g",t_g); results->AddObject(result); + result=new Result(results->Size()+1,time,step,"m_g",m_g); results->AddObject(result); + + //update inputs + model->UpdateInputsFromVector(h_g,ThicknessEnum,VertexEnum); + model->UpdateInputsFromVector(s_g,SurfaceEnum,VertexEnum); + model->UpdateInputsFromVector(b_g,BedEnum,VertexEnum); + + SplitSolutionVectorx(u_g,numberofnodes,3,&vx,&vy,&vz); + model->UpdateInputsFromVector(vx,VxEnum,VertexEnum); + model->UpdateInputsFromVector(vy,VyEnum,VertexEnum); + model->UpdateInputsFromVector(vz,VzEnum,VertexEnum); + xfree((void**)&vx); xfree((void**)&vy); xfree((void**)&vz); + + model->UpdateInputsFromVector(p_g,PressureEnum,VertexEnum); + model->UpdateInputsFromVector(t_g,TemperatureEnum,VertexEnum); + + //Save temporary results every 5 steps + if (step%5==0){ + _printf_("%s"," saving temporary results..."); + + DataSet* processed_results=NULL; + char* outputfilename=NULL; + + model->FindParam(&outputfilename,OutputFileNameEnum); + ProcessResults(&processed_results,results,model,TransientAnalysisEnum); + OutputResults(processed_results,outputfilename); + + delete processed_results; + xfree((void**)&outputfilename); + + _printf_("%s\n"," done."); + } + } + + /*Add analysis_type to results: */ + result=new Result(results->Size()+1,0,1,"analysis_type",EnumAsString(TransientAnalysisEnum)); + results->AddObject(result); + + + /*Free ressources:*/ + VecFree(&u_g); + VecFree(&p_g); + VecFree(&h_g); + VecFree(&s_g); + VecFree(&b_g); + VecFree(&t_g); + VecFree(&m_g); + xfree((void**)&vx); + xfree((void**)&vy); + xfree((void**)&vz); + + //return: + return results; + +}