Index: /issm/trunk/src/c/DataSet/DataSet.cpp =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.cpp (revision 822) +++ /issm/trunk/src/c/DataSet/DataSet.cpp (revision 823) @@ -338,4 +338,33 @@ } +int DataSet::FindResult(Vec* presult,char* name){ + + /*Go through a dataset, and find a Result* object + *whith result name is "name" : */ + + vector::iterator object; + Result* result=NULL; + + int found=0; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + /*Find param type objects: */ + if((*object)->Enum()==ResultEnum()){ + + /*Ok, this object is a result,recover it and ask which name it has: */ + result=(Result*)(*object); + + if (strcmp(result->GetFieldName(),name)==0){ + /*Ok, this is the one! Recover the value of this result: */ + result->GetField(presult); + found=1; + break; + } + } + } + return found; +} + int DataSet::FindResult(void* pvalue, char* name){ Index: /issm/trunk/src/c/DataSet/DataSet.h =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.h (revision 822) +++ /issm/trunk/src/c/DataSet/DataSet.h (revision 823) @@ -44,4 +44,5 @@ int Size(); int FindParam(void* pvalue, char* name); + int FindResult(Vec* presult,char* name); Object* FindParamObject(char* name); void NodeRank(int* ranks); Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 822) +++ /issm/trunk/src/c/Makefile.am (revision 823) @@ -557,4 +557,7 @@ ./parallel/ProcessResults.cpp\ ./parallel/prognostic_core.cpp\ + ./parallel/transient_core.cpp\ + ./parallel/transient_core_2d.cpp\ + ./parallel/transient_core_3d.cpp\ ./parallel/qmu.cpp\ ./parallel/DakotaResponses.cpp\ @@ -567,5 +570,5 @@ bin_PROGRAMS = else -bin_PROGRAMS = diagnostic.exe control.exe thermal.exe prognostic.exe +bin_PROGRAMS = diagnostic.exe control.exe thermal.exe prognostic.exe transient.exe endif @@ -584,2 +587,5 @@ prognostic_exe_CXXFLAGS= -fPIC -D_PARALLEL_ +transient_exe_SOURCES = parallel/transient.cpp +transient_exe_CXXFLAGS= -fPIC -D_PARALLEL_ + Index: /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateParametersPrognostic.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateParametersPrognostic.cpp (revision 822) +++ /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateParametersPrognostic.cpp (revision 823) @@ -26,4 +26,6 @@ double* pressure=NULL; double* thickness=NULL; + double* surface=NULL; + double* bed=NULL; double* accumulation=NULL; double* melting=NULL; @@ -68,4 +70,29 @@ xfree((void**)&thickness); + /*Get surface: */ + ModelFetchData((void**)&surface,NULL,NULL,model_handle,"surface","Matrix","Mat"); + + count++; + param= new Param(count,"s_g",DOUBLEVEC); + if(surface) param->SetDoubleVec(surface,model->numberofnodes,1); + else param->SetDoubleVec(surface,0,0); + parameters->AddObject(param); + + /*Free surface: */ + xfree((void**)&surface); + + /*Get bed: */ + ModelFetchData((void**)&bed,NULL,NULL,model_handle,"bed","Matrix","Mat"); + + count++; + param= new Param(count,"b_g",DOUBLEVEC); + if(bed) param->SetDoubleVec(bed,model->numberofnodes,1); + else param->SetDoubleVec(bed,0,0); + parameters->AddObject(param); + + /*Free bed: */ + xfree((void**)&bed); + + /*Get melting: */ ModelFetchData((void**)&melting,NULL,NULL,model_handle,"melting","Matrix","Mat"); @@ -93,6 +120,4 @@ - - /*Assign output pointer: */ *pparameters=parameters; Index: /issm/trunk/src/c/objects/Result.cpp =================================================================== --- /issm/trunk/src/c/objects/Result.cpp (revision 822) +++ /issm/trunk/src/c/objects/Result.cpp (revision 823) @@ -196,9 +196,12 @@ void Result::GetField(Vec* pfield){ - *pfield=field; + + VecDuplicatePatch(pfield,field); + } void Result::GetField(double** pfield){ - *pfield=dfield; + *pfield=(double*)xmalloc(size*sizeof(double)); + memcpy(*pfield,dfield,size*sizeof(double)); } Index: /issm/trunk/src/c/parallel/SpawnCore.cpp =================================================================== --- /issm/trunk/src/c/parallel/SpawnCore.cpp (revision 822) +++ /issm/trunk/src/c/parallel/SpawnCore.cpp (revision 823) @@ -144,5 +144,5 @@ if(debug)_printf_("Starting transient core\n"); - throw ErrorException(__FUNCT__,"not supported yet!"); + transient_core(results,femmodels,inputs); } Index: /issm/trunk/src/c/parallel/parallel.h =================================================================== --- /issm/trunk/src/c/parallel/parallel.h (revision 822) +++ /issm/trunk/src/c/parallel/parallel.h (revision 823) @@ -19,4 +19,8 @@ void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); void diagnostic_core_linear(Vec* ppg,FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); + +void transient_core(DataSet* results,FemModel* fems, ParameterInputs* inputs); +void transient_core_2d(DataSet* results,FemModel* fems, ParameterInputs* inputs); +void transient_core_3d(DataSet* results,FemModel* fems, ParameterInputs* inputs); //int GradJOrth(WorkspaceParams* workspaceparams); Index: /issm/trunk/src/c/parallel/transient.cpp =================================================================== --- /issm/trunk/src/c/parallel/transient.cpp (revision 823) +++ /issm/trunk/src/c/parallel/transient.cpp (revision 823) @@ -0,0 +1,152 @@ +/*!\file: transient.cpp + * \brief: transient solution + */ + +#include "../issm.h" +#include "./parallel.h" + +#undef __FUNCT__ +#define __FUNCT__ "transient" + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + + +int main(int argc,char* *argv){ + + /*I/O: */ + FILE* fid=NULL; + char* inputfilename=NULL; + char* outputfilename=NULL; + char* lockname=NULL; + char* qmuinname=NULL; + char* qmuoutname=NULL; + char* qmuerrname=NULL; + int numberofnodes; + int qmu_analysis=0; + + /*Fem models : */ + FemModel femmodels[8]; + int dim=-1; + + /*Results: */ + DataSet* results=NULL; + + ParameterInputs* inputs=NULL; + int waitonlock=0; + + /*inputs: */ + double* u_g=NULL; + double* m_g=NULL; + double* a_g=NULL; + double dt; + Param* param=NULL; + + MODULEBOOT(); + + #if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_)) + throw ErrorException(__FUNCT__," parallel executable was compiled without support of parallel libraries!"); + #endif + + PetscInitialize(&argc,&argv,(char *)0,""); + + /*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]; + qmuinname=argv[5]; + qmuoutname=argv[6]; + qmuerrname=argv[7]; + + /*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"); + CreateFemModel(&femmodels[0],fid,"diagnostic","horiz"); + _printf_("\n reading diagnostic vert model data:\n"); + CreateFemModel(&femmodels[1],fid,"diagnostic","vert"); + _printf_("\n reading diagnostic stokes model data:\n"); + CreateFemModel(&femmodels[2],fid,"diagnostic","stokes"); + _printf_("\n reading diagnostic hutter model data:\n"); + CreateFemModel(&femmodels[3],fid,"diagnostic","hutter"); + _printf_("\n reading surface and bed slope computation model data:\n"); + CreateFemModel(&femmodels[4],fid,"slope_compute",""); + _printf_("\n reading prognositc model data:\n"); + CreateFemModel(&femmodels[5],fid,"prognostic",""); + + /*Do we run in 3d?, in which case we need thermal and melting also:*/ + femmodels[0].parameters->FindParam((void*)&dim,"dim"); + if(dim==3){ + _printf_("read and create thermal finite element model:\n"); + CreateFemModel(&femmodels[6],fid,"thermal",NULL); + _printf_("read and create melting finite element model:\n"); + CreateFemModel(&femmodels[7],fid,"melting",NULL); + } + + _printf_("initialize inputs:\n"); + femmodels[5].parameters->FindParam((void*)&numberofnodes,"numberofnodes"); + + femmodels[5].parameters->FindParam((void*)&u_g,"u_g"); + inputs=new ParameterInputs; inputs->Add("velocity",u_g,2,numberofnodes); + + femmodels[5].parameters->FindParam((void*)&m_g,"m_g"); + inputs=new ParameterInputs; inputs->Add("melting",m_g,1,numberofnodes); + + femmodels[5].parameters->FindParam((void*)&a_g,"a_g"); + inputs=new ParameterInputs; inputs->Add("accumulation",a_g,1,numberofnodes); + + femmodels[5].parameters->FindParam((void*)&dt,"dt"); + inputs=new ParameterInputs; inputs->Add("dt",dt); + + _printf_("initialize results:\n"); + results=new DataSet(ResultsEnum()); + + /*are we running the solution sequence, or a qmu wrapper around it? : */ + femmodels[5].parameters->FindParam((void*)&qmu_analysis,"qmu_analysis"); + if(!qmu_analysis){ + + /*run diagnostic analysis: */ + _printf_("call computational core:\n"); + transient_core(results,femmodels,inputs); + + } + else{ + /*run qmu analysis: */ + _printf_("calling qmu analysis on transient core:\n"); + + #ifdef _HAVE_DAKOTA_ + qmu(qmuinname,qmuoutname,qmuerrname,&femmodels[0],inputs,TransientAnalysisEnum(),NoneAnalysisEnum()); + #else + throw ErrorException(__FUNCT__," Dakota not present, cannot do qmu!"); + #endif + } + + _printf_("process results:\n"); + ProcessResults(&results,&femmodels[0],DiagnosticAnalysisEnum()); + + _printf_("write results to disk:\n"); + OutputResults(results,outputfilename); + + _printf_("write lock file:\n"); + femmodels[0].parameters->FindParam((void*)&waitonlock,"waitonlock"); + if (waitonlock){ + WriteLockFile(lockname); + } + + _printf_("closing MPI and Petsc\n"); + PetscFinalize(); + + + /*end module: */ + MODULEEND(); + + return 0; //unix success return; +} Index: /issm/trunk/src/c/parallel/transient_core.cpp =================================================================== --- /issm/trunk/src/c/parallel/transient_core.cpp (revision 823) +++ /issm/trunk/src/c/parallel/transient_core.cpp (revision 823) @@ -0,0 +1,39 @@ +/*!\file: transient_core.cpp + * \brief: core of the transient solution + */ + +#undef __FUNCT__ +#define __FUNCT__ "transient_core" + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./parallel.h" +#include "../issm.h" + +void transient_core(DataSet* results,FemModel* fems, ParameterInputs* inputs){ + + extern int my_rank; + + /*fem models: */ + FemModel* fem_p=NULL; + + int dim=-1; + + + fem_p=fems+5; + + //first recover parameters common to all solutions + fem_p->parameters->FindParam((void*)&dim,"dim"); + + //branch out + if(dim==2){ + transient_core_2d(results,fems,inputs); + } + else if(dim==3){ + transient_core_3d(results,fems,inputs); + } + else throw ErrorException(__FUNCT__,exprintf("%s%i%s"," dimension ",dim," not supported yet!")); + +} Index: /issm/trunk/src/c/parallel/transient_core_2d.cpp =================================================================== --- /issm/trunk/src/c/parallel/transient_core_2d.cpp (revision 823) +++ /issm/trunk/src/c/parallel/transient_core_2d.cpp (revision 823) @@ -0,0 +1,164 @@ +/*!\file: transient_2d_core.cpp + * \brief: core of the transient_2d solution + */ + +#undef __FUNCT__ +#define __FUNCT__ "transient_2d_core" + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./parallel.h" +#include "../issm.h" + +void transient_core_2d(DataSet* results,FemModel* fems, ParameterInputs* inputs){ + + 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: */ + Result* result=NULL; + + /*solutions: */ + int step; + double time; + + Vec u_g=NULL; + double* u_g_serial=NULL; + Vec p_g=NULL; + double* p_g_serial=NULL; + Vec h_g=NULL; + double* h_g_serial=NULL; + Vec s_g=NULL; + double* s_g_serial=NULL; + Vec b_g=NULL; + double* b_g_serial=NULL; + + Vec next_u_g=NULL; + Vec next_p_g=NULL; + Vec next_h_g=NULL; + Vec h_g_intermediary=NULL; + Vec next_s_g=NULL; + Vec next_b_g=NULL; + + /*intermediary: */ + double finaltime; + double dt; + double yts; + DataSet* diagnostic_results=NULL; + DataSet* prognostic_results=NULL; + + + /*flags: */ + int debug=0; + int numberofnodes; + + /*dof recovery: */ + int dof01[2]={0,1}; + int dof2[1]={2}; + int dof012[3]={0,1,2}; + int dof3[1]={3}; + + /*recover fem models: */ + fem_dh=fems+0; + fem_dv=fems+1; + fem_ds=fems+2; + fem_dhu=fems+3; + fem_sl=fems+4; + fem_p=fems+5; + + //first recover parameters common to all solutions + fem_dh->parameters->FindParam((void*)&debug,"debug"); + fem_dh->parameters->FindParam((void*)&finaltime,"ndt"); + fem_dh->parameters->FindParam((void*)&dt,"dt"); + fem_dh->parameters->FindParam((void*)&yts,"yts"); + fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes"); + + /*initialize: */ + step=1; + time=0; + + fem_p->parameters->FindParam((void*)&u_g_serial,"u_g"); + u_g=SerialToVec(u_g_serial,2*numberofnodes); xfree((void**)&u_g_serial); + + fem_p->parameters->FindParam((void*)&p_g_serial,"p_g"); + p_g=SerialToVec(p_g_serial,2*numberofnodes); xfree((void**)&p_g_serial); + + fem_p->parameters->FindParam((void*)&h_g_serial,"h_g"); + h_g=SerialToVec(h_g_serial,2*numberofnodes); xfree((void**)&h_g_serial); + + fem_p->parameters->FindParam((void*)&s_g_serial,"s_g"); + s_g=SerialToVec(s_g_serial,2*numberofnodes); xfree((void**)&s_g_serial); + + fem_p->parameters->FindParam((void*)&b_g_serial,"b_g"); + b_g=SerialToVec(b_g_serial,2*numberofnodes); xfree((void**)&b_g_serial); + + /*Add all these initial vectors in the 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); + + while(timeAdd("thickness",h_g,1,numberofnodes); + inputs->Add("surface",s_g,1,numberofnodes); + inputs->Add("bed",b_g,1,numberofnodes); + + //Deal with velocities. + diagnostic_results=new DataSet(ResultsEnum()); diagnostic_core(diagnostic_results,fems, inputs); + + //extract u_g and p_g from diagnostic_results, and erase diagnostic_results; + diagnostic_results->FindResult(&next_u_g,"u_g"); + diagnostic_results->FindResult(&next_p_g,"p_g"); + delete diagnostic_results; + + //compute new thickness + _printf_("%s\n","computing new thickness"); + + inputs->Add("velocity",next_u_g,2,numberofnodes); + prognostic_results=new DataSet(ResultsEnum()); prognostic_core(prognostic_results,fems, inputs); + + //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(&next_h_g,&next_b_g,&next_s_g, + fem_p->elements, fem_p->nodes,fem_p->loads, fem_p->materials, fem_p->parameters, + h_g_intermediary,b_g,s_g); + VecFree(&h_g_intermediary); + + //plug into results. + result=new Result(results->Size()+1,time,step,"u_g",next_u_g); results->AddObject(result); + result=new Result(results->Size()+1,time,step,"p_g",next_p_g); results->AddObject(result); + result=new Result(results->Size()+1,time,step,"h_g",next_h_g); results->AddObject(result); + result=new Result(results->Size()+1,time,step,"s_g",next_s_g); results->AddObject(result); + result=new Result(results->Size()+1,time,step,"b_g",next_b_g); results->AddObject(result); + + //point our results to next step + VecFree(&u_g); u_g=next_u_g; + VecFree(&p_g); p_g=next_p_g; + VecFree(&h_g); h_g=next_h_g; + VecFree(&s_g); s_g=next_s_g; + VecFree(&b_g); b_g=next_b_g; + + } + +} Index: /issm/trunk/src/c/parallel/transient_core_3d.cpp =================================================================== --- /issm/trunk/src/c/parallel/transient_core_3d.cpp (revision 823) +++ /issm/trunk/src/c/parallel/transient_core_3d.cpp (revision 823) @@ -0,0 +1,20 @@ +/*!\file: transient_3d_core.cpp + * \brief: core of the transient_3d solution + */ + +#undef __FUNCT__ +#define __FUNCT__ "transient_core_3d" + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./parallel.h" +#include "../issm.h" + +void transient_core_3d(DataSet* results,FemModel* fems, ParameterInputs* inputs){ + + throw ErrorException(__FUNCT__," 3d transient not supported yet!"); + + +} Index: /issm/trunk/src/m/solutions/cielo/transient2d.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/transient2d.m (revision 822) +++ /issm/trunk/src/m/solutions/cielo/transient2d.m (revision 823) @@ -33,6 +33,6 @@ solution.p_g=[]; solution.h_g=m_p.parameters.h_g; -solution.s_g=md.surface; -solution.b_g=md.bed; +solution.s_g=m_p.parameters.s_g; +solution.b_g=m_p.parameters.b_g; %initialize inputs