Index: /issm/trunk/src/c/EnumDefinitions/AnalysisTypeAsEnum.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/AnalysisTypeAsEnum.cpp (revision 471) +++ /issm/trunk/src/c/EnumDefinitions/AnalysisTypeAsEnum.cpp (revision 472) @@ -57,4 +57,10 @@ return VertAnalysisEnum(); } + else if (strcmp(analysis_type,"steady")==0){ + return SteadyAnalysisEnum(); + } + else if (strcmp(analysis_type,"transient")==0){ + return TransientAnalysisEnum(); + } else if (strcmp(analysis_type,"")==0){ return NoneAnalysisEnum(); Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp (revision 471) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp (revision 472) @@ -30,4 +30,6 @@ int VertAnalysisEnum(void){ return 33; } int NoneAnalysisEnum(void){ return 34; } +int SteadyAnalysisEnum(void){ return 35; } +int TransientAnalysisEnum(void){ return 36; } /*datasets: */ Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 471) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 472) @@ -49,4 +49,6 @@ int VertAnalysisEnum(void); int NoneAnalysisEnum(void); +int SteadyAnalysisEnum(void); +int TransientAnalysisEnum(void); Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 471) +++ /issm/trunk/src/c/Makefile.am (revision 472) @@ -162,4 +162,6 @@ ./io/FetchRifts.cpp\ ./io/ParameterInputsInit.cpp\ + ./io/pfopen.cpp\ + ./io/pfclose.cpp\ ./EnumDefinitions/EnumDefinitions.h\ ./EnumDefinitions/EnumDefinitions.cpp\ @@ -406,4 +408,6 @@ ./io/FetchRifts.cpp\ ./io/ParameterInputsInit.cpp\ + ./io/pfopen.cpp\ + ./io/pfclose.cpp\ ./EnumDefinitions/EnumDefinitions.h\ ./EnumDefinitions/EnumDefinitions.cpp\ @@ -493,4 +497,5 @@ ./parallel/diagnostic_core_linear.cpp\ ./parallel/diagnostic_core_nonlinear.cpp\ + ./parallel/thermal_core.cpp\ ./parallel/CreateFemModel.cpp\ ./parallel/OutputDiagnostic.cpp\ @@ -498,4 +503,5 @@ ./parallel/control.cpp\ ./parallel/OutputControl.cpp\ + ./parallel/OutputThermal.cpp\ ./parallel/objectivefunctionC.cpp\ ./parallel/GradJCompute.cpp Index: /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp =================================================================== --- /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp (revision 471) +++ /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp (revision 472) @@ -13,5 +13,5 @@ #include "../EnumDefinitions/EnumDefinitions.h" -void PenaltySystemMatricesx(Mat Kgg, Vec pg,DataSet* elements,DataSet* nodes,DataSet* loads,DataSet* materials, +void PenaltySystemMatricesx(Mat Kgg, Vec pg,double* pkmax,DataSet* elements,DataSet* nodes,DataSet* loads,DataSet* materials, int kflag,int pflag,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ @@ -52,5 +52,6 @@ #endif - + /*Assign output pointers:*/ + if(pkmax)*pkmax=kmax; } Index: /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.h =================================================================== --- /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.h (revision 471) +++ /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.h (revision 472) @@ -10,5 +10,5 @@ /* local prototypes: */ -void PenaltySystemMatricesx(Mat Kgg, Vec pg,DataSet* elements,DataSet* nodes,DataSet* loads,DataSet* materials, +void PenaltySystemMatricesx(Mat Kgg, Vec pg,double* pkmax, DataSet* elements,DataSet* nodes,DataSet* loads,DataSet* materials, int kflag,int pflag,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); Index: /issm/trunk/src/c/io/io.h =================================================================== --- /issm/trunk/src/c/io/io.h (revision 471) +++ /issm/trunk/src/c/io/io.h (revision 472) @@ -45,4 +45,8 @@ int FetchRifts(int** priftsnumpenaltypairs,double*** priftspenaltypairs,int** priftsfill,double** priftsfriction,ConstDataHandle model_handle,int numrifts); +/*File I/O: */ +FILE* pfopen(char* filename,char* format); +void* pfclose(FILE* fid,char* filename); + #endif /* _IMDB_H */ Index: /issm/trunk/src/c/io/pfclose.cpp =================================================================== --- /issm/trunk/src/c/io/pfclose.cpp (revision 472) +++ /issm/trunk/src/c/io/pfclose.cpp (revision 472) @@ -0,0 +1,22 @@ +/*!\file: pfclose.cpp + * \brief fclose wrapper for parallel solution + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "stdio.h" +#include "../shared/shared.h" + +#undef __FUNCT__ +#define __FUNCT__ "pfclose" + +void pfclose(FILE* fid,char* filename){ + + /*Close file handle: */ + if(fclose(fid)!=0)throw ErrorException(__FUNCT__,exprintf("%s%s","could not close file ",filename)); +} + Index: /issm/trunk/src/c/io/pfopen.cpp =================================================================== --- /issm/trunk/src/c/io/pfopen.cpp (revision 472) +++ /issm/trunk/src/c/io/pfopen.cpp (revision 472) @@ -0,0 +1,27 @@ +/*!\file: pfopen.cpp + * \brief fopen wrapper for parallel solution + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "stdio.h" +#include "../shared/shared.h" + +#undef __FUNCT__ +#define __FUNCT__ "pfopen" + +FILE* pfopen(char* filename,char* format){ + + FILE* fid=NULL; + + /*Open handle to data on disk: */ + fid=fopen(filename,format); + if(fid==NULL) throw ErrorException(__FUNCT__,exprintf("%s%s%s","could not open file ",filename," for binary reading or writing")); + + return fid; +} + Index: sm/trunk/src/c/parallel/BatchDebug.c =================================================================== --- /issm/trunk/src/c/parallel/BatchDebug.c (revision 471) +++ (revision ) @@ -1,77 +1,0 @@ -/* - * BatchDebug: get a matrix or a vector written to disk, so that it can be read on the matlab - * side and compared with ISSM code: - */ -#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) - -#include "../include/cielo.h" - -#undef __FUNCT__ -#define __FUNCT__ "OutputControl" - - -int BatchDebug(Mat* Kgg,Vec* pg,FemModel* femmodel,char* filename){ - - int noerr=1; - - /* standard output: */ - Vec* tpartition=NULL; - double* serial_tpartition=NULL; - int serial_tpartition_rows; - int dummy=1; - FILE* fid=NULL; - - double* serial_pg=NULL; - int serial_pg_rows; - double* serial_Kgg=NULL; - int serial_Kgg_rows,serial_Kgg_cols; - - /*recover parameters: */ - tpartition=femmodel->tpartition; - - /*serialize partition vector: */ - VecGetSize(*tpartition,&serial_tpartition_rows); - VecToMPISerial(&serial_tpartition,tpartition); - - /*Serialize output matrices and vectors if they exist: */ - if(Kgg){ - MatGetSize(*Kgg,&serial_Kgg_rows,&serial_Kgg_cols); - MatToSerial(&serial_Kgg,Kgg); - } - - if(pg){ - VecGetSize(*pg,&serial_pg_rows); - VecToMPISerial(&serial_pg,pg); - } - - - /* Open output file to write raw binary data: */ - if(my_rank==0){ - fid=fopen(filename,"wb"); - if(fid==NULL){ - _printf_("%s%s\n",__FUNCT__," error message: could not open file ",filename," for binary writing."); - noerr=0; goto cleanup_and_return; - } - - /*Write tpartition: */ - WriteDataToDisk(serial_tpartition,&serial_tpartition_rows,&dummy,"Mat",fid); - - /*Write Kgg if it exists: */ - if(Kgg)WriteDataToDisk(serial_Kgg,&serial_Kgg_rows,&serial_Kgg_cols,"Mat",fid); - - /*Write pg if it exists: */ - if(pg)WriteDataToDisk(serial_pg,&serial_pg_rows,&dummy,"Mat",fid); - - - /*Close file: */ - if(fclose(fid)!=0){ - _printf_("%s%s%s\n",__FUNCT__," error message: could not close file ",filename); - noerr=0; goto cleanup_and_return; - } - } - - cleanup_and_return: - return noerr; -} -#endif //#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) - Index: sm/trunk/src/c/parallel/GradJCheck.c =================================================================== --- /issm/trunk/src/c/parallel/GradJCheck.c (revision 471) +++ (revision ) @@ -1,69 +1,0 @@ -/* - * GradJCheck.c: - */ - -#include "../../../config.h" - -#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) - -#include "../include/cielo.h" -#include "../modules.h" -#include "./parallel.h" - -#undef __FUNCT__ -#define __FUNCT__ "GradJCheck" -#undef CLEANUP -#define CLEANUP GradJCheckLocalCleanup(); - -void GradJCheckLocalCleanup(void); - -int GradJCheck(WorkspaceParams* workspaceparams,int step,int status){ - - int reloop=0; - int i; - - if (step==0){ - /*Ok, this is the first time we are running the GradJSearch. No criterion applied here. : */ - } - else{ - if (workspaceparams->J[step]>workspaceparams->J[step-1]){ - /*Ok, the GradJSearch yielded a worse cost function. reloop and increate maxiter by 5 for the rest of the model. As well, - * lower tolx: */ - if(status==0){ - /*tolerance was breached! Are we under 10^-10?: */ - if (workspaceparams->tolx<10^-10){ - /*Ok, we are never converging!: */ - _printf_("%s\n"," convergence failed. getting out!"); - reloop=0; - } - else{ - workspaceparams->tolx=workspaceparams->tolx/5; - _printf_("%s%g\n"," optimization failed. relooping with tolx changed to: ",workspaceparams->tolx); - reloop=1; - } - } - else{ - /* maxiter was breached. Are we over 50 iterations?: */ - if(workspaceparams->maxiter[step]>50){ - _printf_("%s\n"," convergence failed. getting out!"); - reloop=0; - } - else{ - for(i=step;insteps;i++){ - workspaceparams->maxiter[i]+=5; - } - _printf_("%s%g\n"," optimization failed. relooping with maxiter changed to: ",workspaceparams->maxiter[step]); - reloop=1; - } - } - } - } - return reloop; -} - - -void GradJCheckLocalCleanup(void){ - return; -} -#endif //#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) - Index: sm/trunk/src/c/parallel/GradJOrth.c =================================================================== --- /issm/trunk/src/c/parallel/GradJOrth.c (revision 471) +++ (revision ) @@ -1,57 +1,0 @@ -/* - * GradJOrth.c: - */ - -#include "../../../config.h" - -#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) - -#include "../include/cielo.h" -#include "../modules.h" -#include "./parallel.h" - -#undef __FUNCT__ -#define __FUNCT__ "GradJOrth" -#undef CLEANUP -#define CLEANUP GradJOrthLocalCleanup(); - -void GradJOrthLocalCleanup(void); - -int GradJOrth(WorkspaceParams* workspaceparams){ - - /*Error management: */ - int noerr=1; - int i; - - Vec* gradient=NULL; - Vec* gradient2=NULL; - Vec* oldgradient=NULL; - Vec* newgradient=NULL; - - - for(i=0;inum_control_parameters;i++){ - char* control_type=workspaceparams->control_types[i]; - - gradient=WorkspaceParamsGetParameterGradient(workspaceparams,control_type); - oldgradient=WorkspaceParamsGetOldParameterGradient(workspaceparams,control_type); - - Orthx( &newgradient,gradient,oldgradient); - - /*Save oldgradient: */ - gradient2=xmalloc(sizeof(Vec));VecDuplicate(*gradient,gradient2);VecCopy(*gradient,*gradient2); - WorkspaceParamsSetOldParameterGradient(workspaceparams,gradient2,control_type); - - /*Set new gradient: */ - WorkspaceParamsSetParameterGradient(workspaceparams,newgradient,control_type); - } - - - EXIT(noerr); -} - -void GradJOrthLocalCleanup(void){ - return; -} - -#endif //#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) - Index: sm/trunk/src/c/parallel/GradJSearch.c =================================================================== --- /issm/trunk/src/c/parallel/GradJSearch.c (revision 471) +++ (revision ) @@ -1,328 +1,0 @@ -/* - * GradJSearch.c: - */ - -#include "../../../config.h" - -#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) - -#include "../include/cielo.h" -#include "../modules.h" -#include "./parallel.h" - -#undef __FUNCT__ -#define __FUNCT__ "GradJSearch" -#undef CLEANUP -#define CLEANUP GradJSearchLocalCleanup(); - -void GradJSearchLocalCleanup(void); - - -int GradJSearch(double* search_vector,FemModel* femmodel,int step){ - - /*Error management: */ - int noerr=1; - int i,n; - int dummy; - ParameterInputs* inputs=NULL; - int status; - - //status=GoldenSearch(search_vector,femmodel->workspaceparams->J+step,-1,1,femmodel->workspaceparams->tolx,(int)femmodel->workspaceparams->maxiter[step], - // femmodel->workspaceparams->fit[step],femmodel->workspaceparams->optscal[step],&objectivefunctionC,femmodel); //do only one dimension search for now. - - status=BrentSearch(search_vector,femmodel->workspaceparams->J+step,-1,1,femmodel->workspaceparams->tolx,(int)femmodel->workspaceparams->maxiter[step], - femmodel->workspaceparams->fit[step],femmodel->workspaceparams->optscal[step],&objectivefunctionC,femmodel); //do only one dimension search for now. - - TESTEXIT(noerr); - - return status; -} - -void GradJSearchLocalCleanup(void){ - return; -} - -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*,ParameterInputs*),FemModel* femmodel){ - - double xc, xd, fc, fd; - double oneminustau = 1 - (sqrt(5) - 1) / 2; - int iter = 0; - ParameterInputs* inputs=NULL; - int status; - - inputs=NewParameterInputs(); - - xc = xa + oneminustau * (xb - xa); - fc = (*f)(&xc,fit,optscal,femmodel,inputs); - xd = xb - oneminustau * (xb - xa); - fd = (*f)(&xd,fit,optscal,femmodel,inputs); - do { - iter++; - if (fc < fd) { - xb = xd; - xd = xc; - xc = xa + oneminustau * (xb - xa); - fd = fc; - fc = (*f)(&xc,fit,optscal,femmodel,inputs); - } - else { - xa = xc; - xc = xd; - xd = xb - oneminustau * (xb - xa); - fc = fd; - fd = (*f)(&xd,fit,optscal,femmodel,inputs); - } - _printf_(" iter# %i f(x) %g x %g toler %g/%g iter %i/%i\n",iter,(fc+fd)/2,(xa+xb)/2,fabs(xb-xa),tolerance,iter,maxiter); - } - while (fabs(xb - xa) > tolerance && iter < maxiter); - - if (fabs(xb-xa)tol1){ - - // Yes, so fit parabola - goldenflag=0; - parab_num=(xbest-x1)*(xbest-x1)*(fxbest-fx2)-(xbest-x2)*(xbest-x2)*(fxbest-fx1);; - parab_den=2.0*(xbest-x1)*(fxbest-fx2)-2.0*(xbest-x2)*(fxbest-fx1); - - //reverse p if necessary - if(parab_den>0.0){ - parab_num=-parab_num; - } - parab_den=sqrt(pow(parab_den,2)); - oldintervalgold=intervalgold; - intervalgold=distance; - - // Is the parabola acceptable - if (( sqrt(pow(parab_num,2)) < sqrt(pow(0.5*parab_den*oldintervalgold,2))) && - (parab_num>parab_den*(a-xbest)) && - (parab_num=xm){ - intervalgold=a-xbest; - } - else{ - intervalgold=b-xbest; - } - distance=gold*intervalgold; - } - - // The function must not be evaluated too close to xbest - if(distance<0){ - si=-1; - } - else{ - si=1; - } - if (sqrt(pow(distance,2))>tol1){ - x=xbest+si*sqrt(pow(distance,2)); - } - else{ - x=xbest+si*tol1; - } - - //evaluate function on x - fx = (*f)(&x,fit,optscal,femmodel,inputs); - iter=iter+1; - - // Update a, b, xm, x1, x2, tol1, tol2 - if (fx<=fxbest){ - if (x>=xbest){ - a=xbest; - } - else{ - b=xbest; - } - x1=x2; fx1=fx2; - x2=xbest; fx2=fxbest; - xbest=x; fxbest=fx; - } - - else{ // fx > fxbest - if (x < xbest){ - a=x; - } - else{ - b=x; - } - if ((fx<=fx2) || (x2==xbest)){ - x1=x2; fx1=fx2; - x2=x; fx2=fx; - } - else if ( (fx <= fx1) || (x1 == xbest) || (x1 == x2) ){ - x1=x; fx1=fx; - } - } - xm = 0.5*(a+b); - tol1=seps*pow(pow(xbest,2),0.5)+tolerance/3.0; - tol2=2.0*tol1; - - //print result - if (goldenflag){ - _printf_(" %5i %12.6g %12.6g %12.6g %s\n",iter,x,fx,pow(pow(xbest-xm,2),0.5)," golden"); - } - else{ - _printf_(" %5i %12.6g %12.6g %12.6g %s\n",iter,x,fx,pow(pow(xbest-xm,2),0.5)," parabolic"); - } - - //Stop the optimization? - if (sqrt(pow(xbest-xm,2)) < (tol2-0.5*(b-a))){ - _printf_("\nOptimization terminated:\nthe current x satisfies the termination criteria using 'tolx' of %g \n", tolerance); - loop=0; - status=0; - } - else if (iter>=maxiter){ - _printf_("\nExiting: Maximum number of iterations has been exceeded - increase 'maxiter'\n"); - loop=0; - status=1; - } - else{ - //continue - loop=1; - } - }//end while - - //Now, check that the value on the boundaries are not better than current fxbest - if (fxbest>fxmin){ - xval=xmin; - fval=fxmin; - } - else if (fxbest>fxmax){ - xval=xmax; - fval=fxmax; - } - else{ - xval=xbest; - fval=fxbest; - } - - /*Assign output pointers: */ - *psearch_scalar=xval; - *pJ=fval; - - return status; -} -#endif //#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) - Index: /issm/trunk/src/c/parallel/OutputControl.cpp =================================================================== --- /issm/trunk/src/c/parallel/OutputControl.cpp (revision 471) +++ /issm/trunk/src/c/parallel/OutputControl.cpp (revision 472) @@ -38,6 +38,5 @@ /* Open output file to write raw binary data: */ if(my_rank==0){ - fid=fopen(filename,"wb"); - if(fid==NULL)throw ErrorException(__FUNCT__,exprintf("%s%s%s","could not open file ",filename," for binary writing.")); + fid=pfopen(filename,"wb"); /*Write solution type: */ @@ -62,5 +61,5 @@ /*Close file: */ - if(fclose(fid)!=0)throw ErrorException(__FUNCT__,exprintf("%s%s","could not close file ",filename)); + pfclose(fid,filename); } Index: /issm/trunk/src/c/parallel/OutputDiagnostic.cpp =================================================================== --- /issm/trunk/src/c/parallel/OutputDiagnostic.cpp (revision 471) +++ /issm/trunk/src/c/parallel/OutputDiagnostic.cpp (revision 472) @@ -50,6 +50,5 @@ /* Open output file to write raw binary data: */ if(my_rank==0){ - fid=fopen(filename,"wb"); - if(fid==NULL)throw ErrorException(__FUNCT__,exprintf("%s%s%s","could not open file ",filename," for binary writing.")); + fid=pfopen(filename,"wb"); /*Write solution type: */ @@ -67,5 +66,5 @@ /*Close file: */ - if(fclose(fid)!=0)throw ErrorException(__FUNCT__,exprintf("%s%s","could not close file ",filename)); + pfclose(fid,filename); } Index: sm/trunk/src/c/parallel/OutputThermal.c =================================================================== --- /issm/trunk/src/c/parallel/OutputThermal.c (revision 471) +++ (revision ) @@ -1,81 +1,0 @@ - -/* - OutputThermal.c: output model results for thermal solution. -*/ -#include "../../../config.h" - -#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) - -#include "../include/cielo.h" - -#undef __FUNCT__ -#define __FUNCT__ "OutputThermal" - -#undef CLEANUP -#define CLEANUP OutputThermalLocalCleanup(); - -void OutputThermalLocalCleanup(void); - -int OutputThermal(Vec* t_g,Vec* partition,char* filename,char* analysis_type){ - - /* error handling: */ - int i; - int noerr=1; - - /* output: */ - FILE* fid=NULL; - - /* standard output: */ - double* serial_partition=NULL; - int serial_partition_rows; - int analysis_size; - - double* serial_t_g=NULL; - int serial_t_g_rows; - int dummy=1; - - /*serialize outputs: */ - VecShift(partition,1.0); //matlab indexing - VecGetSize(*partition,&serial_partition_rows); - noerr=VecToMPISerial(&serial_partition,partition);TESTEXIT(noerr); - - VecGetSize(*t_g,&serial_t_g_rows); - noerr=VecToMPISerial(&serial_t_g,t_g);TESTEXIT(noerr); - - /* Open output file to write raw binary data: */ - if(my_rank==0){ - fid=fopen(filename,"wb"); - if(fid==NULL){ - _printf_("%s%s%s%s\n",__FUNCT__," error message: could not open file ",filename," for binary writing."); - noerr=0; goto cleanup_and_return; - } - - /*Write solution type: */ - WriteDataToDisk(analysis_type,NULL,NULL,"String",fid); - - /*Write uset.gsize: */ - WriteDataToDisk(&uset->gsize,NULL,NULL,"Integer",fid); - - /*Write partition: */ - WriteDataToDisk(serial_partition,&serial_partition_rows,&dummy,"Mat",fid); - - /*Write solution to disk: */ - WriteDataToDisk(serial_t_g,&serial_t_g_rows,&dummy,"Mat",fid); - - /*Close file: */ - if(fclose(fid)!=0){ - _printf_("%s%s%s\n",__FUNCT__," error message: could not close file ",filename); - noerr=0; goto cleanup_and_return; - } - } - - cleanup_and_return: - TESTEXIT(noerr); - - EXIT(noerr); -} -void OutputThermalLocalCleanup(void){ - return; -} -#endif //#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) - Index: /issm/trunk/src/c/parallel/OutputThermal.cpp =================================================================== --- /issm/trunk/src/c/parallel/OutputThermal.cpp (revision 472) +++ /issm/trunk/src/c/parallel/OutputThermal.cpp (revision 472) @@ -0,0 +1,104 @@ + +/* + OutputThermal.c: output model results for thermal solution. +*/ +#undef __FUNCT__ +#define __FUNCT__ "OutputThermal" + +#include "../toolkits/toolkits.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../shared/shared.h" +#include "../io/io.h" +#include "../objects/objects.h" + +void OutputThermal(Vec* t_g,Vec* m_g, FemModel* femmodels,char* filename){ + + int i; + extern int my_rank; + + /* output: */ + FILE* fid=NULL; + + /*intermediary: */ + FemModel femmodel; + NodeSets* nodesets=NULL; + Vec partition=NULL; + char* analysis_type="thermal"; + + /* standard output: */ + Vec partition_shifted=NULL; + double* serial_partition=NULL; + + double* serial_tg=NULL; + double* serial_mg=NULL; + + double ndt,dt; + int nsteps; + int sub_analysis_type; + + int one=1; + int gsize; + int nods; + + /*Recover thermal horiz femmodel: */ + femmodel=femmodels[0]; + partition=femmodel.partition; + nodesets=femmodel.nodesets; + + femmodels[0].parameters->FindParam((void*)&dt,"dt"); + femmodels[0].parameters->FindParam((void*)&ndt,"ndt"); + femmodels[0].parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type"); + + /*serialize outputs: */ + VecDuplicate(partition,&partition_shifted); + VecCopy(partition,partition_shifted); + VecShift(partition_shifted,1.0); //matlab indexing starts at 1 + VecToMPISerial(&serial_partition,partition_shifted); + VecGetSize(partition,&nods); + + if(sub_analysis_type==SteadyAnalysisEnum()){ + nsteps=0; + } + else{ + nsteps=(int)(ndt/dt); + } + + /* Open output file to write raw binary data: */ + if(my_rank==0){ + fid=pfopen(filename,"wb"); + + /*Write solution type: */ + WriteDataToDisk(analysis_type,NULL,NULL,"String",fid); + + /*Write uset.gsize: */ + gsize=nodesets->GetGSize(); + WriteDataToDisk(&gsize,NULL,NULL,"Integer",fid); + + /*Write partition: */ + WriteDataToDisk(serial_partition,&nods,&one,"Mat",fid); + + /*Write number of steps: */ + WriteDataToDisk(&nsteps,NULL,NULL,"Integer",fid); + + /*Write solutions to disk: */ + for(i=0;i<=nsteps;i++){ + + xfree((void**)&serial_tg); + xfree((void**)&serial_mg); + VecToMPISerial(&serial_tg,t_g[i]); + VecToMPISerial(&serial_mg,m_g[i]); + + WriteDataToDisk(serial_tg,&gsize,&one,"Mat",fid); + WriteDataToDisk(serial_mg,&gsize,&one,"Mat",fid); + } + + /*Close file: */ + pfclose(fid,filename); + } + + /*Free ressources: */ + VecFree(&partition_shifted); + xfree((void**)&serial_partition); + xfree((void**)&serial_tg); + xfree((void**)&serial_mg); +} Index: sm/trunk/src/c/parallel/ParameterUpdate.c =================================================================== --- /issm/trunk/src/c/parallel/ParameterUpdate.c (revision 471) +++ (revision ) @@ -1,61 +1,0 @@ -/* - * ParameterUpdate.c: - */ - -#include "../../../config.h" - -#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) - -#include "../include/cielo.h" -#include "../modules.h" -#include "./parallel.h" - -#undef __FUNCT__ -#define __FUNCT__ "ParameterUpdate" -#undef CLEANUP -#define CLEANUP ParameterUpdateLocalCleanup(&gradient); - -void ParameterUpdateLocalCleanup(double** pgradient); - -int ParameterUpdate(double* search_vector,int step,WorkspaceParams* workspaceparams,BatchParams* batchparams){ - - /*Error management: */ - int noerr=1; - int i,j; - double* parameter=NULL; - Vec* vec_gradient=NULL; - double* gradient=NULL; - - - //Go through parameters, and update along gradients, multiplying by search_vector. - for(i=0;inum_control_parameters;i++){ - char* control_type=workspaceparams->control_types[i]; - - /*Get parameter, and gradient for the parameter: */ - parameter=WorkspaceParamsGetParameter(workspaceparams,control_type); - vec_gradient=WorkspaceParamsGetParameterGradient(workspaceparams,control_type); - - /*serialize gradient: */ - VecToMPISerial(&gradient,vec_gradient); - - /*Ok, for this parameter, we have a direction. Update the parameter along - * the direction, using the search_vector value: */ - for(j=0;j<(int)(workspaceparams->gsize/6);j++){ - parameter[6*j+0]=parameter[6*j+0]-search_vector[i]*workspaceparams->optscal[step]*gradient[6*j+0]; - } - - /*Now, call on each parameter's self constraining routine: */ - WorkspaceParamsConstrain(workspaceparams,control_type); - } - - /*Done, just return: */ - EXIT(noerr); -} - -void ParameterUpdateLocalCleanup(double** pgradient){ - xfree((void**)pgradient); - return; -} - -#endif //#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) - Index: /issm/trunk/src/c/parallel/cielo/BatchDebug.c =================================================================== --- /issm/trunk/src/c/parallel/cielo/BatchDebug.c (revision 472) +++ /issm/trunk/src/c/parallel/cielo/BatchDebug.c (revision 472) @@ -0,0 +1,77 @@ +/* + * BatchDebug: get a matrix or a vector written to disk, so that it can be read on the matlab + * side and compared with ISSM code: + */ +#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) + +#include "../include/cielo.h" + +#undef __FUNCT__ +#define __FUNCT__ "OutputControl" + + +int BatchDebug(Mat* Kgg,Vec* pg,FemModel* femmodel,char* filename){ + + int noerr=1; + + /* standard output: */ + Vec* tpartition=NULL; + double* serial_tpartition=NULL; + int serial_tpartition_rows; + int dummy=1; + FILE* fid=NULL; + + double* serial_pg=NULL; + int serial_pg_rows; + double* serial_Kgg=NULL; + int serial_Kgg_rows,serial_Kgg_cols; + + /*recover parameters: */ + tpartition=femmodel->tpartition; + + /*serialize partition vector: */ + VecGetSize(*tpartition,&serial_tpartition_rows); + VecToMPISerial(&serial_tpartition,tpartition); + + /*Serialize output matrices and vectors if they exist: */ + if(Kgg){ + MatGetSize(*Kgg,&serial_Kgg_rows,&serial_Kgg_cols); + MatToSerial(&serial_Kgg,Kgg); + } + + if(pg){ + VecGetSize(*pg,&serial_pg_rows); + VecToMPISerial(&serial_pg,pg); + } + + + /* Open output file to write raw binary data: */ + if(my_rank==0){ + fid=fopen(filename,"wb"); + if(fid==NULL){ + _printf_("%s%s\n",__FUNCT__," error message: could not open file ",filename," for binary writing."); + noerr=0; goto cleanup_and_return; + } + + /*Write tpartition: */ + WriteDataToDisk(serial_tpartition,&serial_tpartition_rows,&dummy,"Mat",fid); + + /*Write Kgg if it exists: */ + if(Kgg)WriteDataToDisk(serial_Kgg,&serial_Kgg_rows,&serial_Kgg_cols,"Mat",fid); + + /*Write pg if it exists: */ + if(pg)WriteDataToDisk(serial_pg,&serial_pg_rows,&dummy,"Mat",fid); + + + /*Close file: */ + if(fclose(fid)!=0){ + _printf_("%s%s%s\n",__FUNCT__," error message: could not close file ",filename); + noerr=0; goto cleanup_and_return; + } + } + + cleanup_and_return: + return noerr; +} +#endif //#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) + Index: /issm/trunk/src/c/parallel/cielo/GradJCheck.c =================================================================== --- /issm/trunk/src/c/parallel/cielo/GradJCheck.c (revision 472) +++ /issm/trunk/src/c/parallel/cielo/GradJCheck.c (revision 472) @@ -0,0 +1,69 @@ +/* + * GradJCheck.c: + */ + +#include "../../../config.h" + +#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) + +#include "../include/cielo.h" +#include "../modules.h" +#include "./parallel.h" + +#undef __FUNCT__ +#define __FUNCT__ "GradJCheck" +#undef CLEANUP +#define CLEANUP GradJCheckLocalCleanup(); + +void GradJCheckLocalCleanup(void); + +int GradJCheck(WorkspaceParams* workspaceparams,int step,int status){ + + int reloop=0; + int i; + + if (step==0){ + /*Ok, this is the first time we are running the GradJSearch. No criterion applied here. : */ + } + else{ + if (workspaceparams->J[step]>workspaceparams->J[step-1]){ + /*Ok, the GradJSearch yielded a worse cost function. reloop and increate maxiter by 5 for the rest of the model. As well, + * lower tolx: */ + if(status==0){ + /*tolerance was breached! Are we under 10^-10?: */ + if (workspaceparams->tolx<10^-10){ + /*Ok, we are never converging!: */ + _printf_("%s\n"," convergence failed. getting out!"); + reloop=0; + } + else{ + workspaceparams->tolx=workspaceparams->tolx/5; + _printf_("%s%g\n"," optimization failed. relooping with tolx changed to: ",workspaceparams->tolx); + reloop=1; + } + } + else{ + /* maxiter was breached. Are we over 50 iterations?: */ + if(workspaceparams->maxiter[step]>50){ + _printf_("%s\n"," convergence failed. getting out!"); + reloop=0; + } + else{ + for(i=step;insteps;i++){ + workspaceparams->maxiter[i]+=5; + } + _printf_("%s%g\n"," optimization failed. relooping with maxiter changed to: ",workspaceparams->maxiter[step]); + reloop=1; + } + } + } + } + return reloop; +} + + +void GradJCheckLocalCleanup(void){ + return; +} +#endif //#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) + Index: /issm/trunk/src/c/parallel/cielo/GradJOrth.c =================================================================== --- /issm/trunk/src/c/parallel/cielo/GradJOrth.c (revision 472) +++ /issm/trunk/src/c/parallel/cielo/GradJOrth.c (revision 472) @@ -0,0 +1,57 @@ +/* + * GradJOrth.c: + */ + +#include "../../../config.h" + +#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) + +#include "../include/cielo.h" +#include "../modules.h" +#include "./parallel.h" + +#undef __FUNCT__ +#define __FUNCT__ "GradJOrth" +#undef CLEANUP +#define CLEANUP GradJOrthLocalCleanup(); + +void GradJOrthLocalCleanup(void); + +int GradJOrth(WorkspaceParams* workspaceparams){ + + /*Error management: */ + int noerr=1; + int i; + + Vec* gradient=NULL; + Vec* gradient2=NULL; + Vec* oldgradient=NULL; + Vec* newgradient=NULL; + + + for(i=0;inum_control_parameters;i++){ + char* control_type=workspaceparams->control_types[i]; + + gradient=WorkspaceParamsGetParameterGradient(workspaceparams,control_type); + oldgradient=WorkspaceParamsGetOldParameterGradient(workspaceparams,control_type); + + Orthx( &newgradient,gradient,oldgradient); + + /*Save oldgradient: */ + gradient2=xmalloc(sizeof(Vec));VecDuplicate(*gradient,gradient2);VecCopy(*gradient,*gradient2); + WorkspaceParamsSetOldParameterGradient(workspaceparams,gradient2,control_type); + + /*Set new gradient: */ + WorkspaceParamsSetParameterGradient(workspaceparams,newgradient,control_type); + } + + + EXIT(noerr); +} + +void GradJOrthLocalCleanup(void){ + return; +} + +#endif //#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) + Index: /issm/trunk/src/c/parallel/cielo/GradJSearch.c =================================================================== --- /issm/trunk/src/c/parallel/cielo/GradJSearch.c (revision 472) +++ /issm/trunk/src/c/parallel/cielo/GradJSearch.c (revision 472) @@ -0,0 +1,328 @@ +/* + * GradJSearch.c: + */ + +#include "../../../config.h" + +#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) + +#include "../include/cielo.h" +#include "../modules.h" +#include "./parallel.h" + +#undef __FUNCT__ +#define __FUNCT__ "GradJSearch" +#undef CLEANUP +#define CLEANUP GradJSearchLocalCleanup(); + +void GradJSearchLocalCleanup(void); + + +int GradJSearch(double* search_vector,FemModel* femmodel,int step){ + + /*Error management: */ + int noerr=1; + int i,n; + int dummy; + ParameterInputs* inputs=NULL; + int status; + + //status=GoldenSearch(search_vector,femmodel->workspaceparams->J+step,-1,1,femmodel->workspaceparams->tolx,(int)femmodel->workspaceparams->maxiter[step], + // femmodel->workspaceparams->fit[step],femmodel->workspaceparams->optscal[step],&objectivefunctionC,femmodel); //do only one dimension search for now. + + status=BrentSearch(search_vector,femmodel->workspaceparams->J+step,-1,1,femmodel->workspaceparams->tolx,(int)femmodel->workspaceparams->maxiter[step], + femmodel->workspaceparams->fit[step],femmodel->workspaceparams->optscal[step],&objectivefunctionC,femmodel); //do only one dimension search for now. + + TESTEXIT(noerr); + + return status; +} + +void GradJSearchLocalCleanup(void){ + return; +} + +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*,ParameterInputs*),FemModel* femmodel){ + + double xc, xd, fc, fd; + double oneminustau = 1 - (sqrt(5) - 1) / 2; + int iter = 0; + ParameterInputs* inputs=NULL; + int status; + + inputs=NewParameterInputs(); + + xc = xa + oneminustau * (xb - xa); + fc = (*f)(&xc,fit,optscal,femmodel,inputs); + xd = xb - oneminustau * (xb - xa); + fd = (*f)(&xd,fit,optscal,femmodel,inputs); + do { + iter++; + if (fc < fd) { + xb = xd; + xd = xc; + xc = xa + oneminustau * (xb - xa); + fd = fc; + fc = (*f)(&xc,fit,optscal,femmodel,inputs); + } + else { + xa = xc; + xc = xd; + xd = xb - oneminustau * (xb - xa); + fc = fd; + fd = (*f)(&xd,fit,optscal,femmodel,inputs); + } + _printf_(" iter# %i f(x) %g x %g toler %g/%g iter %i/%i\n",iter,(fc+fd)/2,(xa+xb)/2,fabs(xb-xa),tolerance,iter,maxiter); + } + while (fabs(xb - xa) > tolerance && iter < maxiter); + + if (fabs(xb-xa)tol1){ + + // Yes, so fit parabola + goldenflag=0; + parab_num=(xbest-x1)*(xbest-x1)*(fxbest-fx2)-(xbest-x2)*(xbest-x2)*(fxbest-fx1);; + parab_den=2.0*(xbest-x1)*(fxbest-fx2)-2.0*(xbest-x2)*(fxbest-fx1); + + //reverse p if necessary + if(parab_den>0.0){ + parab_num=-parab_num; + } + parab_den=sqrt(pow(parab_den,2)); + oldintervalgold=intervalgold; + intervalgold=distance; + + // Is the parabola acceptable + if (( sqrt(pow(parab_num,2)) < sqrt(pow(0.5*parab_den*oldintervalgold,2))) && + (parab_num>parab_den*(a-xbest)) && + (parab_num=xm){ + intervalgold=a-xbest; + } + else{ + intervalgold=b-xbest; + } + distance=gold*intervalgold; + } + + // The function must not be evaluated too close to xbest + if(distance<0){ + si=-1; + } + else{ + si=1; + } + if (sqrt(pow(distance,2))>tol1){ + x=xbest+si*sqrt(pow(distance,2)); + } + else{ + x=xbest+si*tol1; + } + + //evaluate function on x + fx = (*f)(&x,fit,optscal,femmodel,inputs); + iter=iter+1; + + // Update a, b, xm, x1, x2, tol1, tol2 + if (fx<=fxbest){ + if (x>=xbest){ + a=xbest; + } + else{ + b=xbest; + } + x1=x2; fx1=fx2; + x2=xbest; fx2=fxbest; + xbest=x; fxbest=fx; + } + + else{ // fx > fxbest + if (x < xbest){ + a=x; + } + else{ + b=x; + } + if ((fx<=fx2) || (x2==xbest)){ + x1=x2; fx1=fx2; + x2=x; fx2=fx; + } + else if ( (fx <= fx1) || (x1 == xbest) || (x1 == x2) ){ + x1=x; fx1=fx; + } + } + xm = 0.5*(a+b); + tol1=seps*pow(pow(xbest,2),0.5)+tolerance/3.0; + tol2=2.0*tol1; + + //print result + if (goldenflag){ + _printf_(" %5i %12.6g %12.6g %12.6g %s\n",iter,x,fx,pow(pow(xbest-xm,2),0.5)," golden"); + } + else{ + _printf_(" %5i %12.6g %12.6g %12.6g %s\n",iter,x,fx,pow(pow(xbest-xm,2),0.5)," parabolic"); + } + + //Stop the optimization? + if (sqrt(pow(xbest-xm,2)) < (tol2-0.5*(b-a))){ + _printf_("\nOptimization terminated:\nthe current x satisfies the termination criteria using 'tolx' of %g \n", tolerance); + loop=0; + status=0; + } + else if (iter>=maxiter){ + _printf_("\nExiting: Maximum number of iterations has been exceeded - increase 'maxiter'\n"); + loop=0; + status=1; + } + else{ + //continue + loop=1; + } + }//end while + + //Now, check that the value on the boundaries are not better than current fxbest + if (fxbest>fxmin){ + xval=xmin; + fval=fxmin; + } + else if (fxbest>fxmax){ + xval=xmax; + fval=fxmax; + } + else{ + xval=xbest; + fval=fxbest; + } + + /*Assign output pointers: */ + *psearch_scalar=xval; + *pJ=fval; + + return status; +} +#endif //#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) + Index: /issm/trunk/src/c/parallel/cielo/ParameterUpdate.c =================================================================== --- /issm/trunk/src/c/parallel/cielo/ParameterUpdate.c (revision 472) +++ /issm/trunk/src/c/parallel/cielo/ParameterUpdate.c (revision 472) @@ -0,0 +1,61 @@ +/* + * ParameterUpdate.c: + */ + +#include "../../../config.h" + +#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) + +#include "../include/cielo.h" +#include "../modules.h" +#include "./parallel.h" + +#undef __FUNCT__ +#define __FUNCT__ "ParameterUpdate" +#undef CLEANUP +#define CLEANUP ParameterUpdateLocalCleanup(&gradient); + +void ParameterUpdateLocalCleanup(double** pgradient); + +int ParameterUpdate(double* search_vector,int step,WorkspaceParams* workspaceparams,BatchParams* batchparams){ + + /*Error management: */ + int noerr=1; + int i,j; + double* parameter=NULL; + Vec* vec_gradient=NULL; + double* gradient=NULL; + + + //Go through parameters, and update along gradients, multiplying by search_vector. + for(i=0;inum_control_parameters;i++){ + char* control_type=workspaceparams->control_types[i]; + + /*Get parameter, and gradient for the parameter: */ + parameter=WorkspaceParamsGetParameter(workspaceparams,control_type); + vec_gradient=WorkspaceParamsGetParameterGradient(workspaceparams,control_type); + + /*serialize gradient: */ + VecToMPISerial(&gradient,vec_gradient); + + /*Ok, for this parameter, we have a direction. Update the parameter along + * the direction, using the search_vector value: */ + for(j=0;j<(int)(workspaceparams->gsize/6);j++){ + parameter[6*j+0]=parameter[6*j+0]-search_vector[i]*workspaceparams->optscal[step]*gradient[6*j+0]; + } + + /*Now, call on each parameter's self constraining routine: */ + WorkspaceParamsConstrain(workspaceparams,control_type); + } + + /*Done, just return: */ + EXIT(noerr); +} + +void ParameterUpdateLocalCleanup(double** pgradient){ + xfree((void**)pgradient); + return; +} + +#endif //#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) + Index: /issm/trunk/src/c/parallel/cielo/thermal_core.c =================================================================== --- /issm/trunk/src/c/parallel/cielo/thermal_core.c (revision 472) +++ /issm/trunk/src/c/parallel/cielo/thermal_core.c (revision 472) @@ -0,0 +1,176 @@ +/* + * cielothermal_core.c: + */ + +#include "../include/cielo.h" +#include "../modules.h" +#include "./parallel.h" + +#undef __FUNCT__ +#define __FUNCT__ "cielothermal_core" + +#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) + +int cielothermal_core(Vec** pt_g,ParameterInputs* inputs,FemModel* femmodel){ + + + + /*Femmodel: */ + BatchParams* batchparams=NULL; + Mat* G_mn=NULL; + Vec* y_s=NULL; + int flag_y_s; + DataSet* bgpdt=NULL; + DataSet* bgpdtb=NULL; + DataSet* est=NULL; + DataSet* lst=NULL; + DataSet* new_lst=NULL; + DataSet* ept=NULL; + DataSet* mpt=NULL; + DataSet* geom3=NULL; + + int noerr=1; + + int dummy; + Vec* t_f=NULL; + Vec* t_g=NULL; + + Mat* K_gg=NULL; + Mat* K_gg_initial=NULL; + Mat* KT_gg=NULL; + Vec* p_g=NULL; + Vec* p_g_initial=NULL; + Mat* K_ff=NULL; + Mat* K_fs=NULL; + Vec* p_f=NULL; + + int kflag,ktflag,pflag; + int converged=0; + int count; + + /*some parameters:*/ + double sparsity; + int connectivity; + int analysis_type_enum; + char* solverstring=NULL; + int debug; + int min_thermal_constraints; + + /*intermediary data: */ + int num_unstable_constraints; + + /*Recover model parameters: */ + batchparams=femmodel->batchparams; + G_mn=femmodel->G_mn; + y_s=femmodel->y_s; + flag_y_s=femmodel->flag_y_s; + bgpdt=femmodel->bgpdt; + bgpdtb=femmodel->bgpdtb; + est=femmodel->est; + lst=femmodel->lst; + ept=femmodel->ept; + mpt=femmodel->mpt; + geom3=femmodel->geom3; + uset=femmodel->uset; //external variable + + sparsity=batchparams->sparsity; + connectivity=batchparams->connectivity; + analysis_type_enum=batchparams->analysis_type_enum; + solverstring=batchparams->solverstring; + debug=batchparams->debug; + min_thermal_constraints=batchparams->min_thermal_constraints; + + + //Start iteration on non-linearity + kflag=1; pflag=1; ktflag=0; //stiffness and load generation only: + converged=0; //flag to break loop on non-linearity + count=1; + + for(;;){ + if(count==1){ + //*Generate system matrices + Emgx(&K_gg, &p_g, &KT_gg,kflag,pflag,ktflag,sparsity,connectivity,bgpdt, bgpdtb, est, lst, ept, mpt, geom3, inputs, analysis_type_enum); + + //Keep initial K_gg, so that penalties don't keep piling up onto K_gg: */ + K_gg_initial=xmalloc(sizeof(Mat)); + MatDuplicate(*K_gg,MAT_COPY_VALUES,K_gg_initial); + + p_g_initial=xmalloc(sizeof(Vec)); + VecDuplicate(*p_g,p_g_initial);VecCopy(*p_g,*p_g_initial); + } + else{ + /*Copy K_gg_initial into K_gg, same for p_g: */ + K_gg=xmalloc(sizeof(Mat)); + MatDuplicate(*K_gg_initial,MAT_COPY_VALUES,K_gg); + + p_g=xmalloc(sizeof(Vec)); + VecDuplicate(*p_g_initial,p_g);VecCopy(*p_g_initial,*p_g); + } + + //Add penalties + PenaltyEmgx(K_gg, p_g, kflag,pflag,bgpdt, bgpdtb, est, lst, ept, mpt, geom3, inputs, analysis_type_enum); + + //Reduce tangent matrix from g size to f size + Reducematrixfromgtof(&K_ff,&dummy,&dummy, &K_fs,&dummy,&dummy, K_gg, G_mn, flag_y_s, + uset->pv_m,uset->msize, uset->pv_n,uset->nsize, uset->pv_f,uset->fsize, uset->pv_s,uset->ssize,uset->gsize,uset->msize,uset->ssize); + + //no need for K_gg and KT_gg anymore + MatFree(&K_gg); MatFree(&KT_gg); + + //Reduce load from g size to f size + Reducerightside(&p_f, p_g, G_mn, K_fs, y_s, flag_y_s, + uset->pv_m, uset->msize, uset->pv_n, uset->nsize, uset->pv_f, uset->fsize, + uset->msize,uset->ssize,uset->fsize); + + //no need for p_g and K_fs anymore + VecFree(&p_g); MatFree(&K_fs); + + //solve + if(uset->fsize>0){ + t_f=xmalloc(sizeof(Vec)); + Solverx(t_f,&dummy,K_ff,dummy,dummy,p_f,dummy,NULL,0,solverstring); + } + else{ + _printf_("All dof are constrained, f_set is empty...\n"); + goto cleanup_and_return; + } + + //no need for K_ff and p_f anymore + MatFree(&K_ff);VecFree(&p_f); + + //Merge back to g set + Mergesolvec( &t_g, t_f, G_mn, y_s, uset->pv_m, uset->msize, uset->pv_n, uset->nsize, uset->pv_f, uset->fsize, uset->pv_s, uset->ssize, + uset->gsize,uset->msize,uset->nsize); + + //t_f not needed anymore + VecFree(&t_f); + + //Deal with penalty loads + ParameterInputsAddFromVec(inputs,t_g,"temperature"); + PenaltyConstraintsx(&new_lst, &converged,&num_unstable_constraints, bgpdt, bgpdtb, est, lst, mpt, inputs,analysis_type_enum); + + //Figure out if convergence is reached. + if(!converged){ + if (debug){ + _printf_("%s %i\n"," #unstable constraints ",num_unstable_constraints); + if (num_unstable_constraints<=min_thermal_constraints){ + converged=1; + } + } + } + + /*Reset lst to new_lst:*/ + DeleteDataSet(&lst); lst=new_lst; + + count++; + if(converged==1)break; + } + + cleanup_and_return: + + /*Assign output pointers: */ + *pt_g=t_g; + return noerr; +} +#endif //#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) + Index: /issm/trunk/src/c/parallel/control.cpp =================================================================== --- /issm/trunk/src/c/parallel/control.cpp (revision 471) +++ /issm/trunk/src/c/parallel/control.cpp (revision 472) @@ -75,6 +75,5 @@ /*Open handle to data on disk: */ - fid=fopen(inputfilename,"rb"); - if(fid==NULL) throw ErrorException(__FUNCT__,exprintf("%s%s%s","could not open file ",inputfilename," for binary reading")); + fid=pfopen(inputfilename,"rb"); _printf_("read and create finite element model:\n"); Index: /issm/trunk/src/c/parallel/diagnostic.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic.cpp (revision 471) +++ /issm/trunk/src/c/parallel/diagnostic.cpp (revision 472) @@ -54,6 +54,5 @@ /*Open handle to data on disk: */ - fid=fopen(inputfilename,"rb"); - if(fid==NULL) throw ErrorException(__FUNCT__,exprintf("%s%s%s","could not open file ",inputfilename," for binary reading")); + fid=pfopen(inputfilename,"rb"); _printf_("read and create finite element model:\n"); Index: /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp (revision 471) +++ /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp (revision 472) @@ -88,5 +88,5 @@ if (debug) _printf_(" Generating penalty matrices\n"); //*Generate penalty system matrices - PenaltySystemMatricesx(Kgg, pg,fem->elements,fem->nodes,loads,fem->materials,kflag,pflag,inputs,analysis_type,sub_analysis_type); + PenaltySystemMatricesx(Kgg, pg,NULL,fem->elements,fem->nodes,loads,fem->materials,kflag,pflag,inputs,analysis_type,sub_analysis_type); if (debug) _printf_(" reducing matrix from g to f set\n"); Index: /issm/trunk/src/c/parallel/parallel.h =================================================================== --- /issm/trunk/src/c/parallel/parallel.h (revision 471) +++ /issm/trunk/src/c/parallel/parallel.h (revision 472) @@ -14,5 +14,5 @@ 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); -int cielothermal_core(Vec** pt_g,ParameterInputs* inputs,FemModel* femmodel); +void thermal_core(Vec* pt_g,double* pmelting_offset,FemModel* femmodel,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); //int GradJOrth(WorkspaceParams* workspaceparams); @@ -29,5 +29,5 @@ //int ParameterUpdate(double* search_vector,int step, WorkspaceParams* workspaceparams,BatchParams* batchparams); void OutputDiagnostic(Vec u_g,Vec p_g, FemModel* femmodels,char* filename); -int OutputThermal(Vec* t_g,Vec* tpartition,char* filename,char* analysis_type,int sub_analysis_type); +void OutputThermal(Vec* t_g,Vec* m_g, FemModel* femmodels,char* filename); void OutputControl(Vec u_g,double* p_g, double* J, int nsteps, Vec partition,char* filename,NodeSets* nodesets); void WriteLockFile(char* filename); Index: /issm/trunk/src/c/parallel/thermal.cpp =================================================================== --- /issm/trunk/src/c/parallel/thermal.cpp (revision 471) +++ /issm/trunk/src/c/parallel/thermal.cpp (revision 472) @@ -17,4 +17,6 @@ int main(int argc,char* *argv){ + int i,n; + /*I/O: */ FILE* fid=NULL; @@ -30,12 +32,19 @@ Vec u_g=NULL; Vec p_g=NULL; + double dt; + double ndt; + int nsteps; + int debug=0; - /*solution vectors: */ - Vec t_g=NULL; - Vec m_g=NULL; + Vec* t_g=NULL; + Vec* m_g=NULL; ParameterInputs* inputs=NULL; + Param* param=NULL; + int waitonlock=0; + int sub_analysis_type; + double melting_offset; MODULEBOOT(); @@ -51,5 +60,4 @@ MPI_Comm_size(MPI_COMM_WORLD,&num_procs); - _printf_("recover , input file name and output file name:\n"); inputfilename=argv[2]; outputfilename=argv[3]; @@ -57,9 +65,9 @@ /*Open handle to data on disk: */ - fid=fopen(inputfilename,"rb"); - if(fid==NULL) throw ErrorException(__FUNCT__,exprintf("%s%s%s","could not open file ",inputfilename," for binary reading")); + fid=pfopen(inputfilename,"rb"); - _printf_("read and create finite element model:\n"); + _printf_("read and create thermal finite element model:\n"); CreateFemModel(&femmodels[0],fid,"thermal",""); + _printf_("read and create melting finite element model:\n"); CreateFemModel(&femmodels[1],fid,"melting",""); @@ -68,17 +76,69 @@ femmodels[0].parameters->FindParam((void*)&p_g,"p_g"); femmodels[0].parameters->FindParam((void*)&numberofnodes,"numberofnodes"); + femmodels[0].parameters->FindParam((void*)&dt,"dt"); + femmodels[0].parameters->FindParam((void*)&ndt,"ndt"); + femmodels[0].parameters->FindParam((void*)&waitonlock,"waitonlock"); + femmodels[0].parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type"); + femmodels[0].parameters->FindParam((void*)&debug,"debug"); inputs=new ParameterInputs; - //inputs->Add("velocity",u_g_initial,3,numberofnodes); + inputs->Add("velocity",u_g,3,numberofnodes); + inputs->Add("pressure",p_g,1,numberofnodes); + inputs->Add("dt",dt); - //erase velocities: - //param=(Param*)femmodels[0].parameters->FindParamObject("u_g"); - //femmodels[0].parameters->DeleteObject((Object*)param); + //erase velocities and pressure embedded in parameters + param=(Param*)femmodels[0].parameters->FindParamObject("u_g"); + femmodels[0].parameters->DeleteObject((Object*)param); + param=(Param*)femmodels[0].parameters->FindParamObject("p_g"); + femmodels[0].parameters->DeleteObject((Object*)param); + param=(Param*)femmodels[1].parameters->FindParamObject("u_g"); + femmodels[1].parameters->DeleteObject((Object*)param); + param=(Param*)femmodels[1].parameters->FindParamObject("p_g"); + femmodels[1].parameters->DeleteObject((Object*)param); - _printf_("call computational core:\n"); - diagnostic_core(&u_g,&p_g,femmodels,inputs); + if(sub_analysis_type==SteadyAnalysisEnum()){ + if(debug)_printf_("computing temperatures:\n"); + thermal_core(&t_g[0],&melting_offset,&femmodels[0],inputs,ThermalAnalysisEnum(),SteadyAnalysisEnum()); + + inputs->Add("temperature",t_g[0],1,numberofnodes); + inputs->Add("melting_offset",melting_offset); + + if(debug)_printf_("computing melting:\n"); + diagnostic_core_linear(&m_g[0],&femmodels[1],inputs,MeltingAnalysisEnum(),SteadyAnalysisEnum()); + } + else{ + + nsteps=(int)(ndt/dt); + + /*allocate t_g and m_g arrays: */ + t_g=(Vec*)xmalloc((nsteps+1)*sizeof(Vec)); + m_g=(Vec*)xmalloc((nsteps+1)*sizeof(Vec)); + + //initialize temperature and melting + femmodels[0].parameters->FindParam((void*)&t_g[0],"t_g"); + femmodels[1].parameters->FindParam((void*)&m_g[0],"m_g"); + + //erase temperature and melting embedded in parameters + param=(Param*)femmodels[0].parameters->FindParamObject("t_g"); + femmodels[0].parameters->DeleteObject((Object*)param); + param=(Param*)femmodels[1].parameters->FindParamObject("m_g"); + femmodels[1].parameters->DeleteObject((Object*)param); + + for(i=0;iAdd("temperature",t_g[i],1,numberofnodes); + thermal_core(&t_g[i+1],&melting_offset,&femmodels[0],inputs,ThermalAnalysisEnum(),TransientAnalysisEnum()); + + if(debug)_printf_("computing melting:\n"); + inputs->Add("temperature",t_g[i+1],1,numberofnodes); + inputs->Add("melting_offset",melting_offset); + diagnostic_core_linear(&m_g[0],&femmodels[1],inputs,MeltingAnalysisEnum(),TransientAnalysisEnum()); + } + } _printf_("write results to disk:\n"); - OutputDiagnostic(u_g,p_g,&femmodels[0],outputfilename); + OutputThermal(t_g,m_g,&femmodels[0],outputfilename); _printf_("write lock file:\n"); Index: sm/trunk/src/c/parallel/thermal_core.c =================================================================== --- /issm/trunk/src/c/parallel/thermal_core.c (revision 471) +++ (revision ) @@ -1,176 +1,0 @@ -/* - * cielothermal_core.c: - */ - -#include "../include/cielo.h" -#include "../modules.h" -#include "./parallel.h" - -#undef __FUNCT__ -#define __FUNCT__ "cielothermal_core" - -#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) - -int cielothermal_core(Vec** pt_g,ParameterInputs* inputs,FemModel* femmodel){ - - - - /*Femmodel: */ - BatchParams* batchparams=NULL; - Mat* G_mn=NULL; - Vec* y_s=NULL; - int flag_y_s; - DataSet* bgpdt=NULL; - DataSet* bgpdtb=NULL; - DataSet* est=NULL; - DataSet* lst=NULL; - DataSet* new_lst=NULL; - DataSet* ept=NULL; - DataSet* mpt=NULL; - DataSet* geom3=NULL; - - int noerr=1; - - int dummy; - Vec* t_f=NULL; - Vec* t_g=NULL; - - Mat* K_gg=NULL; - Mat* K_gg_initial=NULL; - Mat* KT_gg=NULL; - Vec* p_g=NULL; - Vec* p_g_initial=NULL; - Mat* K_ff=NULL; - Mat* K_fs=NULL; - Vec* p_f=NULL; - - int kflag,ktflag,pflag; - int converged=0; - int count; - - /*some parameters:*/ - double sparsity; - int connectivity; - int analysis_type_enum; - char* solverstring=NULL; - int debug; - int min_thermal_constraints; - - /*intermediary data: */ - int num_unstable_constraints; - - /*Recover model parameters: */ - batchparams=femmodel->batchparams; - G_mn=femmodel->G_mn; - y_s=femmodel->y_s; - flag_y_s=femmodel->flag_y_s; - bgpdt=femmodel->bgpdt; - bgpdtb=femmodel->bgpdtb; - est=femmodel->est; - lst=femmodel->lst; - ept=femmodel->ept; - mpt=femmodel->mpt; - geom3=femmodel->geom3; - uset=femmodel->uset; //external variable - - sparsity=batchparams->sparsity; - connectivity=batchparams->connectivity; - analysis_type_enum=batchparams->analysis_type_enum; - solverstring=batchparams->solverstring; - debug=batchparams->debug; - min_thermal_constraints=batchparams->min_thermal_constraints; - - - //Start iteration on non-linearity - kflag=1; pflag=1; ktflag=0; //stiffness and load generation only: - converged=0; //flag to break loop on non-linearity - count=1; - - for(;;){ - if(count==1){ - //*Generate system matrices - Emgx(&K_gg, &p_g, &KT_gg,kflag,pflag,ktflag,sparsity,connectivity,bgpdt, bgpdtb, est, lst, ept, mpt, geom3, inputs, analysis_type_enum); - - //Keep initial K_gg, so that penalties don't keep piling up onto K_gg: */ - K_gg_initial=xmalloc(sizeof(Mat)); - MatDuplicate(*K_gg,MAT_COPY_VALUES,K_gg_initial); - - p_g_initial=xmalloc(sizeof(Vec)); - VecDuplicate(*p_g,p_g_initial);VecCopy(*p_g,*p_g_initial); - } - else{ - /*Copy K_gg_initial into K_gg, same for p_g: */ - K_gg=xmalloc(sizeof(Mat)); - MatDuplicate(*K_gg_initial,MAT_COPY_VALUES,K_gg); - - p_g=xmalloc(sizeof(Vec)); - VecDuplicate(*p_g_initial,p_g);VecCopy(*p_g_initial,*p_g); - } - - //Add penalties - PenaltyEmgx(K_gg, p_g, kflag,pflag,bgpdt, bgpdtb, est, lst, ept, mpt, geom3, inputs, analysis_type_enum); - - //Reduce tangent matrix from g size to f size - Reducematrixfromgtof(&K_ff,&dummy,&dummy, &K_fs,&dummy,&dummy, K_gg, G_mn, flag_y_s, - uset->pv_m,uset->msize, uset->pv_n,uset->nsize, uset->pv_f,uset->fsize, uset->pv_s,uset->ssize,uset->gsize,uset->msize,uset->ssize); - - //no need for K_gg and KT_gg anymore - MatFree(&K_gg); MatFree(&KT_gg); - - //Reduce load from g size to f size - Reducerightside(&p_f, p_g, G_mn, K_fs, y_s, flag_y_s, - uset->pv_m, uset->msize, uset->pv_n, uset->nsize, uset->pv_f, uset->fsize, - uset->msize,uset->ssize,uset->fsize); - - //no need for p_g and K_fs anymore - VecFree(&p_g); MatFree(&K_fs); - - //solve - if(uset->fsize>0){ - t_f=xmalloc(sizeof(Vec)); - Solverx(t_f,&dummy,K_ff,dummy,dummy,p_f,dummy,NULL,0,solverstring); - } - else{ - _printf_("All dof are constrained, f_set is empty...\n"); - goto cleanup_and_return; - } - - //no need for K_ff and p_f anymore - MatFree(&K_ff);VecFree(&p_f); - - //Merge back to g set - Mergesolvec( &t_g, t_f, G_mn, y_s, uset->pv_m, uset->msize, uset->pv_n, uset->nsize, uset->pv_f, uset->fsize, uset->pv_s, uset->ssize, - uset->gsize,uset->msize,uset->nsize); - - //t_f not needed anymore - VecFree(&t_f); - - //Deal with penalty loads - ParameterInputsAddFromVec(inputs,t_g,"temperature"); - PenaltyConstraintsx(&new_lst, &converged,&num_unstable_constraints, bgpdt, bgpdtb, est, lst, mpt, inputs,analysis_type_enum); - - //Figure out if convergence is reached. - if(!converged){ - if (debug){ - _printf_("%s %i\n"," #unstable constraints ",num_unstable_constraints); - if (num_unstable_constraints<=min_thermal_constraints){ - converged=1; - } - } - } - - /*Reset lst to new_lst:*/ - DeleteDataSet(&lst); lst=new_lst; - - count++; - if(converged==1)break; - } - - cleanup_and_return: - - /*Assign output pointers: */ - *pt_g=t_g; - return noerr; -} -#endif //#if defined(_PARALLEL_) && defined(_HAVE_PETSC_) - Index: /issm/trunk/src/c/parallel/thermal_core.cpp =================================================================== --- /issm/trunk/src/c/parallel/thermal_core.cpp (revision 472) +++ /issm/trunk/src/c/parallel/thermal_core.cpp (revision 472) @@ -0,0 +1,109 @@ +/*!\file: thermal_core.cpp + * \brief: core of the thermal solution + */ + +#undef __FUNCT__ +#define __FUNCT__ "thermal_core" + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../issm.h" + +void thermal_core(Vec* ptg,double* pmelting_offset,FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ + + /*solution : */ + Vec tg=NULL; + Vec tf=NULL; + double melting_offset; + + /*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_thermal_constraints; + + /*parameters:*/ + int kflag,pflag,connectivity,numberofdofspernode; + char* solver_string=NULL; + int debug=0; + + /*Recover parameters: */ + kflag=1; pflag=1; + + fem->parameters->FindParam((void*)&connectivity,"connectivity"); + fem->parameters->FindParam((void*)&numberofdofspernode,"numberofdofspernode"); + fem->parameters->FindParam((void*)&numberofnodes,"numberofnodes"); + fem->parameters->FindParam((void*)&solver_string,"solverstring"); + fem->parameters->FindParam((void*)&debug,"debug"); + fem->parameters->FindParam((void*)&min_thermal_constraints,"min_thermal_constraints"); + + count=1; + converged=0; + for(;;){ + + if(debug)_printf_("%s\n","starting direct shooting method"); + + /*Update parameters: */ + UpdateFromInputsx(fem->elements,fem->nodes,fem->loads, fem->materials,inputs); + + //*Generate system matrices + SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type,sub_analysis_type); + //apply penalties + PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,inputs,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 (debug) _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(debug)_printf_("%s\n","solving"); + Solverx(&tf, Kff, pf, tf, solver_string); + + //no need for Kff and pf anymore + MatFree(&Kff);VecFree(&pf); + + if (debug) _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 (debug) _printf_(" penalty constraints\n"); + inputs->Add("temperature",tg,numberofdofspernode,numberofnodes); + + PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,fem->loads,fem->materials,inputs,analysis_type,sub_analysis_type); + + if (!converged){ + if(debug)_printf_("%s%i\n","#unstable constraints",num_unstable_constraints); + + if (num_unstable_constraints