Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 8428) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 8429) @@ -49,8 +49,5 @@ HydrologyAnalysisEnum, //transient - Transient2DSolutionEnum, - Transient2DAnalysisEnum, - Transient3DSolutionEnum, - Transient3DAnalysisEnum, + TransientSolutionEnum, GroundingLineMigration2DSolutionEnum, SteadyAnalysisEnum, Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 8428) +++ /issm/trunk/src/c/Makefile.am (revision 8429) @@ -1312,6 +1312,5 @@ ./solutions/hydrology_core.cpp\ ./solutions/hydrology_core_step.cpp\ - ./solutions/transient2d_core.cpp\ - ./solutions/transient3d_core.cpp\ + ./solutions/transient_core.cpp\ ./solutions/groundinglinemigration2d_core.cpp\ ./solutions/steadystate_core.cpp\ Index: /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp =================================================================== --- /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 8428) +++ /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 8429) @@ -44,8 +44,5 @@ case HydrologySolutionEnum : return "HydrologySolution"; case HydrologyAnalysisEnum : return "HydrologyAnalysis"; - case Transient2DSolutionEnum : return "Transient2DSolution"; - case Transient2DAnalysisEnum : return "Transient2DAnalysis"; - case Transient3DSolutionEnum : return "Transient3DSolution"; - case Transient3DAnalysisEnum : return "Transient3DAnalysis"; + case TransientSolutionEnum : return "TransientSolution"; case GroundingLineMigration2DSolutionEnum : return "GroundingLineMigration2DSolution"; case SteadyAnalysisEnum : return "SteadyAnalysis"; Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 8428) +++ /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 8429) @@ -18,8 +18,8 @@ void CreateDataSets(Elements** pelements,Nodes** pnodes, Vertices** pvertices, Materials** pmaterials, Constraints** pconstraints, Loads** ploads,Parameters** pparameters,IoModel* iomodel,FILE* iomodel_handle,const int solution_type,const int analysis_type,const int nummodels,int analysis_counter){ - bool continuous = true; - Elements *elements = NULL; - Materials *materials = NULL; - Parameters *parameters = NULL; + bool continuous = true; + Elements *elements = NULL; + Materials *materials = NULL; + Parameters *parameters = NULL; /*Create elements, vertices and materials, independent of analysis_type: */ @@ -32,5 +32,5 @@ /*Now, branch onto analysis dependent model generation: */ switch(analysis_type){ - case DiagnosticHorizAnalysisEnum: + case DiagnosticHorizAnalysisEnum: CreateNodesDiagnosticHoriz(pnodes, iomodel,iomodel_handle); CreateConstraintsDiagnosticHoriz(pconstraints,iomodel,iomodel_handle); @@ -39,5 +39,5 @@ break; - case DiagnosticVertAnalysisEnum: + case DiagnosticVertAnalysisEnum: CreateNodesDiagnosticVert(pnodes, iomodel,iomodel_handle); CreateConstraintsDiagnosticVert(pconstraints,iomodel,iomodel_handle); @@ -120,5 +120,5 @@ /*Update Elements in case we are running a transient solution: */ parameters=*pparameters; - if(analysis_counter==(nummodels-1)&& (solution_type==Transient2DSolutionEnum || solution_type==Transient3DSolutionEnum)){ + if(analysis_counter==(nummodels-1)&& solution_type==TransientSolutionEnum){ UpdateElementsTransient(elements,parameters,iomodel,iomodel_handle,analysis_counter,analysis_type); } @@ -130,4 +130,3 @@ * will need to start at the end of the updated counters: */ UpdateCounters(iomodel,pnodes,ploads,pconstraints); - } Index: /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 8428) +++ /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 8429) @@ -22,16 +22,14 @@ /*output: */ - Elements* elements=NULL; - Nodes* nodes=NULL; - Vertices* vertices=NULL; - Materials* materials=NULL; - Constraints* constraints=NULL; - Loads* loads=NULL; - Parameters* parameters=NULL; + Elements *elements = NULL; + Nodes *nodes = NULL; + Vertices *vertices = NULL; + Materials *materials = NULL; + Constraints *constraints = NULL; + Loads *loads = NULL; + Parameters *parameters = NULL; - /*intermediary: */ - IoModel* iomodel=NULL; - - iomodel = new IoModel(IOMODEL); + /*Initialize IoModel from input file*/ + IoModel* iomodel = new IoModel(IOMODEL); SetVerbosityLevel(iomodel->verbose); @@ -39,4 +37,8 @@ analysis_type=analysis_type_list[i]; + + /*Hack for trasient runs (to be improved)*/ + if(solution_type==TransientSolutionEnum && analysis_type==ThermalAnalysisEnum && iomodel->dim==2) continue; + if(solution_type==TransientSolutionEnum && analysis_type==MeltingAnalysisEnum && iomodel->dim==2) continue; _printf_(VerboseMProcessor()," create datasets for analysis %s\n",EnumToStringx(analysis_type)); Index: /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp =================================================================== --- /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 8428) +++ /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 8429) @@ -42,8 +42,5 @@ else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum; else if (strcmp(name,"HydrologyAnalysis")==0) return HydrologyAnalysisEnum; - else if (strcmp(name,"Transient2DSolution")==0) return Transient2DSolutionEnum; - else if (strcmp(name,"Transient2DAnalysis")==0) return Transient2DAnalysisEnum; - else if (strcmp(name,"Transient3DSolution")==0) return Transient3DSolutionEnum; - else if (strcmp(name,"Transient3DAnalysis")==0) return Transient3DAnalysisEnum; + else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum; else if (strcmp(name,"GroundingLineMigration2DSolution")==0) return GroundingLineMigration2DSolutionEnum; else if (strcmp(name,"SteadyAnalysis")==0) return SteadyAnalysisEnum; Index: /issm/trunk/src/c/shared/Matlab/PrintfFunction.cpp =================================================================== --- /issm/trunk/src/c/shared/Matlab/PrintfFunction.cpp (revision 8428) +++ /issm/trunk/src/c/shared/Matlab/PrintfFunction.cpp (revision 8429) @@ -6,4 +6,6 @@ #include "stdarg.h" #include "stdio.h" +#include "../../shared/shared.h" +#include "../../include/include.h" #ifdef _SERIAL_ @@ -12,30 +14,51 @@ int PrintfFunction(char* format,...){ + /*http://linux.die.net/man/3/vsnprintf*/ /*string to be printed: */ - char string[1024]; //assume that nobody will print more than 1024 characters at once. - int string_size; - extern int my_rank; - extern int num_procs; + char *buffer = NULL; + int n,size = 100; + int string_size; + extern int my_rank; + extern int num_procs; //variable list of arguments - va_list ap; + va_list args; - //First use vsprintf to get the whole input string. - va_start(ap,format); - #ifndef WIN32 - string_size=vsprintf(string,format,ap); //printf style coding - #else - string_size=vsprintf_s(string,format,ap); //printf style coding - #endif - va_end(ap); + while(true){ - //Ok, if we are running in paralle, get node 0 to print - #if defined(_PARALLEL_) - if(my_rank==0)printf(string); - #else - mexPrintf(string); - #endif + /*allocate buffer for given string size*/ + buffer=(char*)xmalloc(size*sizeof(char)); + /* Try to print in the allocated space. */ + va_start(args, format); +#ifndef WIN32 + n=vsnprintf(buffer,size,format,args); +#else + n=vsnprintf_s(buffer,size,format,args); +#endif + va_end(args); + + /* If that worked, return the string. */ + if(n>-1 && n-1) /* glibc 2.1 */ + size=n+1; /* precisely what is needed */ + else /* glibc 2.0 */ + size*=2; /* twice the old size */ + + xfree((void**)&buffer); + } + + /*Ok, if we are running in paralle, get node 0 to print*/ +#if defined(_PARALLEL_) + if(my_rank==0)printf(buffer); +#else + mexPrintf(buffer); +#endif + + /*Clean up and return*/ + xfree((void**)&buffer); return 1; } Index: /issm/trunk/src/c/solutions/CorePointerFromSolutionEnum.cpp =================================================================== --- /issm/trunk/src/c/solutions/CorePointerFromSolutionEnum.cpp (revision 8428) +++ /issm/trunk/src/c/solutions/CorePointerFromSolutionEnum.cpp (revision 8429) @@ -49,9 +49,6 @@ solutioncore=&bedslope_core; break; - case Transient2DSolutionEnum: - solutioncore=&transient2d_core; - break; - case Transient3DSolutionEnum: - solutioncore=&transient3d_core; + case TransientSolutionEnum: + solutioncore=&transient_core; break; case GroundingLineMigration2DSolutionEnum: Index: /issm/trunk/src/c/solutions/SolutionConfiguration.cpp =================================================================== --- /issm/trunk/src/c/solutions/SolutionConfiguration.cpp (revision 8428) +++ /issm/trunk/src/c/solutions/SolutionConfiguration.cpp (revision 8429) @@ -17,5 +17,4 @@ #include "../include/include.h" #include "../solvers/solvers.h" - void SolutionConfiguration(int** panalyses,int* pnumanalyses, void (**psolutioncore)(FemModel*),int solutiontype){ @@ -99,16 +98,5 @@ break; - case Transient2DSolutionEnum: - numanalyses=6; - analyses=(int*)xmalloc(numanalyses*sizeof(int)); - analyses[0]=DiagnosticHorizAnalysisEnum; - analyses[1]=DiagnosticVertAnalysisEnum; - analyses[2]=DiagnosticHutterAnalysisEnum; - analyses[3]=SurfaceSlopeAnalysisEnum; - analyses[4]=BedSlopeAnalysisEnum; - analyses[5]=PrognosticAnalysisEnum; - break; - - case Transient3DSolutionEnum: + case TransientSolutionEnum: numanalyses=8; analyses=(int*)xmalloc(numanalyses*sizeof(int)); @@ -118,7 +106,7 @@ analyses[3]=SurfaceSlopeAnalysisEnum; analyses[4]=BedSlopeAnalysisEnum; - analyses[5]=PrognosticAnalysisEnum; - analyses[6]=ThermalAnalysisEnum; - analyses[7]=MeltingAnalysisEnum; + analyses[5]=ThermalAnalysisEnum; + analyses[6]=MeltingAnalysisEnum; + analyses[7]=PrognosticAnalysisEnum; break; @@ -140,4 +128,3 @@ else xfree((void**)&analyses); if(psolutioncore)*psolutioncore=solutioncore; - } Index: /issm/trunk/src/c/solutions/controltao_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/controltao_core.cpp (revision 8428) +++ /issm/trunk/src/c/solutions/controltao_core.cpp (revision 8429) @@ -17,38 +17,49 @@ /*TAO*/ int i,n,info; - TaoMethod method = "tao_blmvm"; - //TaoMethod method = "tao_lmvm"; - //TaoMethod method = "tao_cg"; - //TaoMethod method = "tao_gpcg"; -> Hessian + //TaoMethod method = "tao_blmvm"; + //TaoMethod method = "tao_lmvm"; + TaoMethod method = "tao_cg"; TaoTerminateReason reason; TAO_SOLVER tao; TAO_APPLICATION controlapp; - Vec initial_solution=NULL; - AppCtx user; /* user-defined work context */ + Vec initial_solution = NULL; + AppCtx user; + PetscInt iter; + double ff,gnorm; -#ifdef _HAVE_TAO_ + /*Initialize TAO*/ int argc; char **args; PetscGetArgs(&argc,&args); int ierr=TaoInitialize(&argc,&args,(char*)0,""); if(ierr) _error_("Could not initialize Tao"); -#endif + + /*Line search options*/ + info = PetscOptionsSetValue("-tao_ls_stepmax","10e11"); if(info) _error_("STOP"); //does not work + info = PetscOptionsSetValue("-tao_ls_stepmin","10e5"); if(info) _error_("STOP"); //does not work + info = PetscOptionsSetValue("-tao_ls_maxfev","8"); if(info) _error_("STOP"); + + /*TAO options: http://www.mcs.anl.gov/research/projects/tao/docs/manualpages/solver/TaoSetFromOptions.html*/ + info = PetscOptionsSetValue("-tao_monitor",""); if(info) _error_("STOP"); + info = PetscOptionsSetValue("-tao_gatol","10e-18"); if(info) _error_("STOP"); + info = PetscOptionsSetValue("-tao_max_its","10"); if(info) _error_("STOP"); + info = PetscOptionsSetValue("-tao_max_funcs","20"); if(info) _error_("STOP"); + + /*Additional options*/ + //info = PetscOptionsSetValue("-info","/u/astrid-r1b/morlighe/svn/issm/trunk/test/NightlyRun/taolog.txt"); if(info) _error_("STOP"); /*Initialize argument*/ user.femmodel=femmodel; + /*Set up and solve TAO*/ + GetVectorFromInputsx(&initial_solution,femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,RheologyBbarEnum,VertexEnum); info = TaoCreate(PETSC_COMM_WORLD,method,&tao); if(info) _error_("STOP"); info = TaoApplicationCreate(PETSC_COMM_WORLD,&controlapp); if(info) _error_("STOP"); - GetVectorFromInputsx(&initial_solution,femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,RheologyBbarEnum,VertexEnum); info = TaoAppSetInitialSolutionVec(controlapp,initial_solution); if(info) _error_("STOP"); info = TaoAppSetObjectiveAndGradientRoutine(controlapp,FormFunctionGradient,(void*)&user); if(info) _error_("STOP"); info = TaoSetOptions(controlapp,tao); if(info) _error_("STOP"); - info = TaoSetTolerances(tao,1e-24,1e-24,1e-24,1e-24); if(info) _error_("STOP"); - //info = TaoSetFunctionLowerBound(tao,0.77); if(info) _error_("STOP"); info = TaoSolveApplication(controlapp,tao); //if(info) _error_("STOP"); - PetscInt iter; - double ff,gnorm; + + /*Get solution status*/ info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,0,0,&reason); //CHKERRQ(info); - - switch(reason){ - /*http://www.mcs.anl.gov/research/projects/tao/docs/manualpages/solver/TaoGetTerminationReason.html*/ + switch(reason){ /*http://www.mcs.anl.gov/research/projects/tao/docs/manualpages/solver/TaoGetTerminationReason.html*/ case TAO_CONVERGED_ATOL: _printf_(true,"TAO_CONVERGED_ATOL (res <= atol)\n"); break; case TAO_CONVERGED_RTOL: _printf_(true,"TAO_CONVERGED_RTOL (res/res0 <= rtol)\n"); break; @@ -72,8 +83,8 @@ } info = TaoView(tao); if(info) _error_("STOP"); + + /*Clean up*/ info = TaoDestroy(tao); if(info) _error_("STOP"); info = TaoAppDestroy(controlapp); if(info) _error_("STOP"); - - /*Clean up*/ VecFree(&initial_solution); @@ -89,13 +100,25 @@ Vec gradient = NULL; + /*Temp*/ +// double* Xserial=NULL; +// VecToMPISerial(&Xserial,X); +// printf("X= [%20.20g %20.20g %20.20g]\n",Xserial[0],Xserial[1],Xserial[2]); + /*End Temp*/ + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,SurfaceAbsVelMisfitEnum,CmResponseEnum); InputUpdateFromVectorx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,X,RheologyBbarEnum,VertexEnum); adjointdiagnostic_core(user->femmodel); Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, RheologyBbarEnum); - //VecView(gradient,PETSC_VIEWER_STDOUT_SELF); - //VecScale(gradient,-1.0); + VecScale(gradient,10e7); + //VecScale(gradient,-1.); VecCopy(gradient,G); + CostFunctionx(fcn, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters,SurfaceAbsVelMisfitEnum); + + //printf("X\n"); + //VecView(X,PETSC_VIEWER_STDOUT_SELF); + + //printf("Gradient\n"); //VecView(G,PETSC_VIEWER_STDOUT_SELF); - CostFunctionx(fcn, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters,SurfaceAbsVelMisfitEnum); + printf("f(x) = %g\n",*fcn); return 0; Index: /issm/trunk/src/c/solutions/issm.cpp =================================================================== --- /issm/trunk/src/c/solutions/issm.cpp (revision 8428) +++ /issm/trunk/src/c/solutions/issm.cpp (revision 8429) @@ -12,13 +12,13 @@ /*I/O: */ - FILE *output_fid = NULL; - FILE *petscoptionsfid = NULL; - char *lockname = NULL; - bool qmu_analysis = false; - bool control_analysis = false; - bool waitonlock = false; + FILE *output_fid = NULL; + FILE *petscoptionsfid = NULL; + char *lockname = NULL; + bool qmu_analysis = false; + bool control_analysis = false; + bool waitonlock = false; /*FemModel: */ - FemModel *femmodel = NULL; + FemModel *femmodel = NULL; /*configuration: */ Index: /issm/trunk/src/c/solutions/solutions.h =================================================================== --- /issm/trunk/src/c/solutions/solutions.h (revision 8428) +++ /issm/trunk/src/c/solutions/solutions.h (revision 8429) @@ -30,6 +30,5 @@ void slopecompute_core(FemModel* femmodel); void steadystate_core(FemModel* femmodel); -void transient2d_core(FemModel* femmodel); -void transient3d_core(FemModel* femmodel); +void transient_core(FemModel* femmodel); void groundinglinemigration2d_core(FemModel* femmodel); double objectivefunctionC(double search_scalar,OptArgs* optargs); Index: sm/trunk/src/c/solutions/transient2d_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/transient2d_core.cpp (revision 8428) +++ (revision ) @@ -1,98 +1,0 @@ -/*!\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 "./solutions.h" -#include "../modules/modules.h" -#include - -void transient2d_core(FemModel* femmodel){ - - int i; - - /*parameters: */ - double finaltime; - double dt ,yts; - int dim = -1; - int solution_type; - bool control_analysis,isdiagnostic,isprognostic,isthermal; - bool time_adapt=false; - int output_frequency; - int gl_migration; - - /*intermediary: */ - int step; - double time; - - /* recover parameters: */ - femmodel->parameters->FindParam(&dim,DimEnum); - femmodel->parameters->FindParam(&finaltime,NdtEnum); - femmodel->parameters->FindParam(&dt,DtEnum); - femmodel->parameters->FindParam(&yts,YtsEnum); - femmodel->parameters->FindParam(&control_analysis,ControlAnalysisEnum); - femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); - femmodel->parameters->FindParam(&output_frequency,OutputFrequencyEnum); - femmodel->parameters->FindParam(&time_adapt,TimeAdaptEnum); - femmodel->parameters->FindParam(&gl_migration,GroundingLineMigrationEnum); - femmodel->parameters->FindParam(&isdiagnostic,IsDiagnosticEnum); - femmodel->parameters->FindParam(&isprognostic,IsPrognosticEnum); - femmodel->parameters->FindParam(&isthermal,IsThermalEnum); - - /*initialize: */ - step=0; - time=0; - - while(time < finaltime - (yts*DBL_EPSILON)){ //make sure we run up to finaltime. - - /*Increment*/ - if(time_adapt){ - TimeAdaptx(&dt,femmodel->elements, femmodel->nodes,femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); - if(time+dt>finaltime) dt=finaltime-time; - femmodel->parameters->SetParam(dt,DtEnum); - } - time+=dt; - femmodel->parameters->SetParam(time,TimeEnum); - step+=1; - - _printf_(VerboseSolution(),"%s%g%s%i%s%g%s%g\n","time [yr]: ",time/yts," iteration number: ",step,"/",floor((finaltime-time)/dt)," dt [yr]: ",dt/yts); - - if(isdiagnostic){ - _printf_(VerboseSolution(),"%s\n"," computing new velocity"); - diagnostic_core(femmodel); - } - - if(isprognostic){ - _printf_(VerboseSolution(),"%s\n"," computing new thickness"); - prognostic_core(femmodel); - } - - if (gl_migration!=NoneEnum){ - _printf_(VerboseSolution(),"%s\n"," computing new grounding line position"); - GroundingLineMigrationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); - } - - if(solution_type==Transient2DSolutionEnum && !control_analysis && (step%output_frequency==0 || time==finaltime)){ - _printf_(VerboseSolution(),"%s\n"," saving results\n"); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VelEnum,step,time); - if(dim==3) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceMassBalanceEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BasalMeltingRateEnum,step,time); - if(gl_migration!=NoneEnum)InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ElementOnIceShelfEnum,step,time); - - /*unload results*/ - _printf_(VerboseSolution(),"%s"," saving temporary results"); - OutputResultsx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,&femmodel->results,step,time); - } - } - -} Index: sm/trunk/src/c/solutions/transient3d_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/transient3d_core.cpp (revision 8428) +++ (revision ) @@ -1,95 +1,0 @@ -/*!\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 "./solutions.h" -#include "../modules/modules.h" -#include "../include/include.h" -#include "../solvers/solvers.h" -#include - -void transient3d_core(FemModel* femmodel){ - - /*parameters: */ - double finaltime; - double dt,yts; - bool control_analysis,isdiagnostic,isprognostic,isthermal; - bool time_adapt=false; - int solution_type; - int output_frequency; - - /*intermediary: */ - int step; - double time; - - //first recover parameters common to all solutions - femmodel->parameters->FindParam(&finaltime,NdtEnum); - femmodel->parameters->FindParam(&dt,DtEnum); - femmodel->parameters->FindParam(&yts,YtsEnum); - femmodel->parameters->FindParam(&control_analysis,ControlAnalysisEnum); - femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); - femmodel->parameters->FindParam(&output_frequency,OutputFrequencyEnum); - femmodel->parameters->FindParam(&time_adapt,TimeAdaptEnum); - femmodel->parameters->FindParam(&isdiagnostic,IsDiagnosticEnum); - femmodel->parameters->FindParam(&isprognostic,IsPrognosticEnum); - femmodel->parameters->FindParam(&isthermal,IsThermalEnum); - - /*initialize: */ - step=0; - time=0; - - while(time < finaltime - (yts*DBL_EPSILON)){ //make sure we run up to finaltime. - - /*Increment*/ - if(time_adapt){ - TimeAdaptx(&dt,femmodel->elements, femmodel->nodes,femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); - if(time+dt>finaltime) dt=finaltime-time; - femmodel->parameters->SetParam(dt,DtEnum); - } - step+=1; - time+=dt; - femmodel->parameters->SetParam(time,TimeEnum); - - _printf_(VerboseSolution(),"%s%g%s%i%s%g%s%g\n","time [yr]: ",time/yts," iteration number: ",step,"/",floor((finaltime-time)/dt)," dt [yr]: ",dt/yts); - - if(isthermal){ - _printf_(VerboseSolution()," computing temperatures:\n"); - thermal_core_step(femmodel,step,time); - } - - if(isdiagnostic){ - _printf_(VerboseSolution(),"%s\n"," computing new velocity"); - diagnostic_core(femmodel); - } - - if(isprognostic){ - _printf_(VerboseSolution(),"%s\n"," computing new thickness"); - prognostic_core(femmodel); - _printf_(VerboseSolution(),"%s\n"," updating vertices positions"); - UpdateVertexPositionsx(femmodel->elements, femmodel->nodes,femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); - } - - if(solution_type==Transient3DSolutionEnum && !control_analysis && (step%output_frequency==0 || time==finaltime)){ - _printf_(VerboseSolution(),"%s\n"," saving results\n"); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VelEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BasalMeltingRateEnum,step,time); - - /*unload results*/ - _printf_(VerboseSolution(),"%s"," saving temporary results"); - OutputResultsx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,&femmodel->results,step,time); - } - } - -} Index: /issm/trunk/src/c/solutions/transient_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/transient_core.cpp (revision 8429) +++ /issm/trunk/src/c/solutions/transient_core.cpp (revision 8429) @@ -0,0 +1,102 @@ +/*!\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 "./solutions.h" +#include "../modules/modules.h" +#include "../include/include.h" +#include "../solvers/solvers.h" +#include + +void transient_core(FemModel* femmodel){ + + /*parameters: */ + double finaltime,dt,yts; + bool control_analysis,isdiagnostic,isprognostic,isthermal; + bool time_adapt=false; + int solution_type; + int output_frequency; + int dim,gl_migration; + + /*intermediary: */ + int step; + double time; + + //first recover parameters common to all solutions + femmodel->parameters->FindParam(&dim,DimEnum); + femmodel->parameters->FindParam(&finaltime,NdtEnum); + femmodel->parameters->FindParam(&dt,DtEnum); + femmodel->parameters->FindParam(&yts,YtsEnum); + femmodel->parameters->FindParam(&control_analysis,ControlAnalysisEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + femmodel->parameters->FindParam(&output_frequency,OutputFrequencyEnum); + femmodel->parameters->FindParam(&time_adapt,TimeAdaptEnum); + femmodel->parameters->FindParam(&gl_migration,GroundingLineMigrationEnum); + femmodel->parameters->FindParam(&isdiagnostic,IsDiagnosticEnum); + femmodel->parameters->FindParam(&isprognostic,IsPrognosticEnum); + femmodel->parameters->FindParam(&isthermal,IsThermalEnum); + + /*initialize: */ + step=0; + time=0; + + while(time < finaltime - (yts*DBL_EPSILON)){ //make sure we run up to finaltime. + + /*Increment*/ + if(time_adapt){ + TimeAdaptx(&dt,femmodel->elements, femmodel->nodes,femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); + if(time+dt>finaltime) dt=finaltime-time; + femmodel->parameters->SetParam(dt,DtEnum); + } + step+=1; + time+=dt; + femmodel->parameters->SetParam(time,TimeEnum); + + _printf_(VerboseSolution(),"%s%g%s%i%s%g%s%g\n","time [yr]: ",time/yts," iteration number: ",step,"/",floor((finaltime-time)/dt)," dt [yr]: ",dt/yts); + + if(isthermal && dim==3){ + _printf_(VerboseSolution()," computing temperatures:\n"); + thermal_core_step(femmodel,step,time); + } + + if(isdiagnostic){ + _printf_(VerboseSolution(),"%s\n"," computing new velocity"); + diagnostic_core(femmodel); + } + + if(isprognostic){ + _printf_(VerboseSolution(),"%s\n"," computing new thickness"); + prognostic_core(femmodel); + _printf_(VerboseSolution(),"%s\n"," updating vertices positions"); + UpdateVertexPositionsx(femmodel->elements, femmodel->nodes,femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); + } + + if (gl_migration!=NoneEnum){ + if(dim==3) _error_("Grounding line migration not implemented in 3d"); + _printf_(VerboseSolution(),"%s\n"," computing new grounding line position"); + GroundingLineMigrationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + } + + if(solution_type==TransientSolutionEnum && !control_analysis && (step%output_frequency==0 || time==finaltime)){ + _printf_(VerboseSolution(),"%s\n"," saving results\n"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,step,time); + if(dim==3) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VelEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time); + if(dim==3) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,step,time); + if(dim==3) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BasalMeltingRateEnum,step,time); + + /*unload results*/ + _printf_(VerboseSolution(),"%s"," saving temporary results"); + OutputResultsx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,&femmodel->results,step,time); + } + } +}