Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 4028) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 4029) @@ -20,7 +20,11 @@ VerticesEnum, /*}}}*/ + /*Solution types {{{1 */ + SolutionTypeEnum, + DiagnosticSolutionEnum, + SurfaceSlopeSolutionEnum, + BedSlopeSolutionEnum, + /*}}}*/ /*Analysis types {{{1 */ - SolutionTypeEnum, - AnalysisEnum, AnalysisTypeEnum, SubAnalysisTypeEnum, @@ -46,11 +50,9 @@ SteadyAnalysisEnum, //slope - SurfaceAnalysisEnum, SlopeComputeAnalysisEnum, - BedXAnalysisEnum, - BedYAnalysisEnum, - BedAnalysisEnum, - SurfaceXAnalysisEnum, - SurfaceYAnalysisEnum, + BedSlopeXAnalysisEnum, + BedSlopeYAnalysisEnum, + SurfaceSlopeXAnalysisEnum, + SurfaceSlopeYAnalysisEnum, //prognostic Balancedthickness2AnalysisEnum, Index: /issm/trunk/src/c/solutions/diagnostic_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/diagnostic_core.cpp (revision 4028) +++ /issm/trunk/src/c/solutions/diagnostic_core.cpp (revision 4029) @@ -22,5 +22,5 @@ double stokesreconditioning; - //first recover parameters needed to drive the solution + /* recover parameters: */ femmodel->parameters->FindParam(&verbose,VerboseEnum); femmodel->parameters->FindParam(&dim,DimEnum); @@ -40,17 +40,11 @@ /*Compute slopes: */ - if(ishutter){ - slope_core(femmodel,SurfaceXAnalysisEnum); - slope_core(femmodel,SurfaceYAnalysisEnum); - } - if(isstokes){ - slope_core(femmodel,BedSlopeXAnalysisEnum); - slope_core(femmodel,BedSlopeYAnalysisEnum); - } + if(ishutter) surfaceslope_core(femmodel); + if(isstokes) bedslope_core(femmodel); if(ishutter){ if(verbose)_printf_("%s\n"," computing hutter velocities..."); - solve_linear(NULL,femmodel,DiagnosticAnalysisEnum,HutterAnalysisEnum); + solver_linear(NULL,femmodel,DiagnosticAnalysisEnum,HutterAnalysisEnum); if (ismacayealpattyn) ResetBoundaryConditions(femmodel,DiagnosticAnalysisEnum,HorizAnalysisEnum); @@ -60,5 +54,5 @@ if(verbose)_printf_("%s\n"," computing horizontal velocities..."); - diagnostic_core_nonlinear(NULL,NULL,NULL,femmodel,false,DiagnosticAnalysisEnum,HorizAnalysisEnum); //false means we expose loads to changes inside the solution + solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,false,DiagnosticAnalysisEnum,HorizAnalysisEnum); //false means we expose loads to changes inside the solution } @@ -67,5 +61,5 @@ if(verbose)_printf_("%s\n"," computing vertical velocities..."); - solve_linear(NULL,femmodel,DiagnosticAnalysisEnum,VertAnalysisEnum); + solver_linear(NULL,femmodel,DiagnosticAnalysisEnum,VertAnalysisEnum); if (isstokes){ @@ -79,5 +73,5 @@ if(verbose)_printf_("%s\n"," computing stokes velocities and pressure ..."); - diagnostic_core_nonlinear(NULL,NULL,NULL,NULL,fem_ds,DiagnosticAnalysisEnum,StokesAnalysisEnum); + solver_diagnostic_nonlinear(NULL,NULL,NULL,NULL,fem_ds,DiagnosticAnalysisEnum,StokesAnalysisEnum); } } Index: /issm/trunk/src/c/solutions/solutions.h =================================================================== --- /issm/trunk/src/c/solutions/solutions.h (revision 4028) +++ /issm/trunk/src/c/solutions/solutions.h (revision 4029) @@ -3,6 +3,6 @@ */ -#ifndef PARALLEL_H_ -#define PARALLEL_H_ +#ifndef SOLUTIONS_H_ +#define SOLUTIONS_H_ #include "../objects/objects.h" @@ -29,11 +29,4 @@ void slope_core(FemModel* fem,int sub_analysis_type); -/*computational cores: */ -void thermal_core_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem,int analysis_type,int sub_analysis_type); -void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* fem,bool conserve_loads,int analysis_type,int sub_analysis_type); -void diagnostic_core_linear(Vec* ppg,FemModel* fem,int analysis_type,int sub_analysis_type); -void solver_linear(Vec* pug, FemModel* femmodel,int analysis_type,int sub_analysis_type); - - //int GradJOrth(WorkspaceParams* workspaceparams); void convergence(int* pconverged, Mat K_ff,Vec p_f,Vec u_f,Vec u_f_old,Parameters* parameters); Index: /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp (revision 4029) +++ /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp (revision 4029) @@ -0,0 +1,170 @@ +/*!\file: solver_diagnostic_nonlinear.cpp + * \brief: core of the diagnostic solution for non linear materials + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../modules/modules.h" +#include "./solutions.h" + +void solver_diagnostic_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, FemModel* fem,bool conserve_loads, int analysis_type,int sub_analysis_type){ + + + /*solution : */ + Vec ug=NULL; + Vec uf=NULL; + Vec old_ug=NULL; + Vec old_uf=NULL; + DataSet* loads=NULL; + + /*intermediary: */ + Mat Kgg=NULL; + Mat Kff=NULL; + Mat Kfs=NULL; + Vec pg=NULL; + Vec pf=NULL; + int converged; + int constraints_converged; + int num_unstable_constraints; + int count; + int numberofnodes; + int min_mechanical_constraints; + int max_nonlinear_iterations; + + /*parameters:*/ + int kflag,pflag; + char* solver_string=NULL; + int verbose=0; + + //Set active analysis_type: + fem->SetActiveAnalysis(analysis_type); + + /*Recover parameters: */ + kflag=1; pflag=1; + fem->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + fem->parameters->FindParam(&solver_string,SolverStringEnum); + fem->parameters->FindParam(&verbose,VerboseEnum); + fem->parameters->FindParam(&min_mechanical_constraints,MinMechanicalConstraintsEnum); + fem->parameters->FindParam(&max_nonlinear_iterations,MaxNonlinearIterationsEnum); + + /*Were loads requested as output? : */ + if(conserve_loads){ + loads=fem->loads->Copy(); //protect loads from being modified by the solution + } + else{ + loads=fem->loads; //modify loads in this solution + } + + count=1; + converged=0; + + /*Start non-linear iteration using input velocity: */ + GetSolutionFromInputsx(&ug, fem->elements, fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters, analysis_type, sub_analysis_type); + Reducevectorgtofx(&uf, ug, fem->nodesets); + + for(;;){ + + //save pointer to old velocity + VecFree(&old_ug);old_ug=ug; + VecFree(&old_uf);old_uf=uf; + + if (verbose) _printf_(" Generating matrices\n"); + //*Generate system matrices + SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); + + if (verbose) _printf_(" Generating penalty matrices\n"); + //*Generate penalty system matrices + PenaltySystemMatricesx(Kgg, pg,NULL,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); + + if (verbose) _printf_(" reducing matrix from g to f set\n"); + /*!Reduce matrix from g to f size:*/ + Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets); + + /*Free ressources: */ + MatFree(&Kgg); + + if (verbose) _printf_(" reducing load from g to f set\n"); + /*!Reduce load from g to f size: */ + Reduceloadfromgtofx(&pf, pg, fem->Gmn, Kfs, fem->ys, fem->nodesets); + + //no need for pg and Kfs anymore + VecFree(&pg); + MatFree(&Kfs); + + /*Solve: */ + if (verbose) _printf_(" solving\n"); + Solverx(&uf, Kff, pf, old_uf, solver_string); + + //Merge back to g set + if (verbose) _printf_(" merging solution from f to g set\n"); + Mergesolutionfromftogx(&ug, uf,fem->Gmn,fem->ys,fem->nodesets); + + //Update inputs using new solution: + UpdateInputsFromSolutionx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,ug,analysis_type, sub_analysis_type); + + //Deal with penalty loads + if (verbose) _printf_(" penalty constraints\n"); + PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,analysis_type,sub_analysis_type); + + if(verbose)_printf_(" number of unstable constraints: %i\n",num_unstable_constraints); + + /*Figure out if convergence is reached.*/ + convergence(&converged,Kff,pf,uf,old_uf,fem->parameters); + MatFree(&Kff);VecFree(&pf); + + /*add converged to inputs: */ + UpdateInputsFromConstantx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,converged,ConvergedEnum); + + //rift convergence + if (!constraints_converged) { + if (converged){ + if (num_unstable_constraints <= min_mechanical_constraints) converged=1; + else converged=0; + } + } + + /*Increase count: */ + count++; + if(converged==1)break; + if(count>=max_nonlinear_iterations){ + _printf_(" maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations); + break; + } + } + + /*extrude if we are in 3D: */ + if (dim==3){ + if(verbose)_printf_("%s\n","extruding velocity and pressure in 3d..."); + InputsExtrudex( fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,VxEnum,0); + InputsExtrudex( fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,VyEnum,0); + } + + //more output might be needed, when running in control + if(pKff0){ + + kflag=1; pflag=0; //stiffness generation only + + SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); + Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets); + MatFree(&Kgg);VecFree(&pg); + + } + + /*Delete loads only if no ouput was requested: */ + if(!conserve_loads)delete loads; + + /*clean up*/ + VecFree(&uf); + VecFree(&old_uf); + VecFree(&old_ug); + xfree((void**)&solver_string); + + /*Assign output pointers: */ + if(pug)*pug=ug; + else VecFree(&ug); + + if(pKff0)*pKff0=Kff; + if(pKfs0)*pKfs0=Kfs; + +} Index: /issm/trunk/src/c/solvers/solver_linear.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_linear.cpp (revision 4029) +++ /issm/trunk/src/c/solvers/solver_linear.cpp (revision 4029) @@ -0,0 +1,65 @@ +/*!\file: solver_linear.cpp + * \brief: numerical core of linear solutions + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../modules/modules.h" + +void solver_linear(Vec* pug, FemModel* fem,int analysis_type,int sub_analysis_type){ + + /*parameters:*/ + int kflag,pflag; + int verbose=0; + char* solver_string=NULL; + + /*output: */ + Vec ug=NULL; + Vec uf=NULL; + + /*intermediary: */ + Mat Kgg=NULL; + Mat Kff=NULL; + Mat Kfs=NULL; + Vec pg=NULL; + Vec pf=NULL; + + /*Recover parameters: */ + kflag=1; pflag=1; + fem->parameters->FindParam(&verbose,VerboseEnum); + fem->parameters->FindParam(&solver_string,SolverStringEnum); + + //*Generate system matrices + if (verbose) _printf_(" Generating matrices\n"); + SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); + + if (verbose) _printf_(" Generating penalty matrices\n"); + //*Generate penalty system matrices + PenaltySystemMatricesx(Kgg, pg,NULL,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); + + /*!Reduce matrix from g to f size:*/ + if (verbose) _printf_(" reducing matrix from g to f set\n"); + Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets); MatFree(&Kgg); + + /*!Reduce load from g to f size: */ + if (verbose) _printf_(" reducing load from g to f set\n"); + Reduceloadfromgtofx(&pf, pg, fem->Gmn, Kfs, fem->ys, fem->nodesets);VecFree(&pg); MatFree(&Kfs); + + /*Solve: */ + if (verbose) _printf_(" solving\n"); + Solverx(&uf, Kff, pf, NULL, solver_string); MatFree(&Kff); VecFree(&pf); + + //Merge back to g set + if (verbose) _printf_(" merging solution from f to g set\n"); + Mergesolutionfromftogx(&ug, uf,fem->Gmn,fem->ys,fem->nodesets);VecFree(&uf); + + //Update inputs using new solution: + UpdateInputsFromSolutionx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,ug,analysis_type, sub_analysis_type); + + /*free ressources: */ + xfree((void**)&solver_string); + + /*Assign output pointers:*/ + if(pug)*pug=ug; +} Index: /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp (revision 4029) +++ /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp (revision 4029) @@ -0,0 +1,135 @@ +/*!\file: solver_thermal_nonlinear.cpp + * \brief: core of the thermal solution + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../modules/modules.h" + +void solver_thermal_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem,int analysis_type,int sub_analysis_type){ + + /*solution : */ + Vec tg=NULL; + Vec tf=NULL; + Vec tf_old=NULL; + double melting_offset; + + /*intermediary: */ + Mat Kgg=NULL; + Mat Kgg_nopenalty=NULL; + Mat Kff=NULL; + Mat Kfs=NULL; + Vec pg=NULL; + Vec pg_nopenalty=NULL; + Vec pf=NULL; + + int converged; + int constraints_converged; + int num_unstable_constraints; + int count; + int numberofnodes; + int min_thermal_constraints; + bool reset_penalties; + + /*parameters:*/ + int kflag,pflag; + char* solver_string=NULL; + int verbose=0; + bool lowmem=0; + + /*Recover parameters: */ + kflag=1; pflag=1; + + fem->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + fem->parameters->FindParam(&solver_string,SolverStringEnum); + fem->parameters->FindParam(&verbose,VerboseEnum); + fem->parameters->FindParam(&lowmem,LowmemEnum); + fem->parameters->FindParam(&min_thermal_constraints,MinThermalConstraintsEnum); + + count=1; + converged=0; + + for(;;){ + + if(verbose)_printf_("%s\n","starting direct shooting method"); + + if(count==1) reset_penalties=1; else reset_penalties=0; + UpdateInputsFromConstantx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,reset_penalties,ResetPenaltiesEnum); + + //*Generate system matrices + if (!lowmem){ + + /*Compute Kgg_nopenalty and pg_nopenalty once for all: */ + if (count==1){ + SystemMatricesx(&Kgg_nopenalty, &pg_nopenalty,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); + } + + /*Copy K_gg_nopenalty into Kgg, same for pg: */ + MatDuplicate(Kgg_nopenalty,MAT_COPY_VALUES,&Kgg); + VecDuplicatePatch(&pg,pg_nopenalty); + + //apply penalties each time + PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); + } + else{ + SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); + //apply penalties + PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); + } + + /*!Reduce matrix from g to f size:*/ + Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets); + + /*Free ressources: */ + MatFree(&Kgg); + + if (verbose) _printf_(" reducing load from g to f set\n"); + /*!Reduce load from g to f size: */ + Reduceloadfromgtofx(&pf, pg, fem->Gmn, Kfs, fem->ys, fem->nodesets); + + //no need for pg and Kfs anymore + VecFree(&pg); + MatFree(&Kfs); + + /*Solve: */ + if(verbose)_printf_("%s\n","solving"); + VecFree(&tf); + Solverx(&tf, Kff, pf,tf_old, solver_string); + VecFree(&tf_old); VecDuplicatePatch(&tf_old,tf); + + //no need for Kff and pf anymore + MatFree(&Kff);VecFree(&pf);VecFree(&tg); + + if (verbose) _printf_(" merging solution from f to g set\n"); + //Merge back to g set + Mergesolutionfromftogx(&tg, tf,fem->Gmn,fem->ys,fem->nodesets); + + //Deal with penalty loads + if (verbose) _printf_(" penalty constraints\n"); + PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,analysis_type,sub_analysis_type); + + //Update inputs using new solution: + UpdateInputsFromVectorx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,tg,TemperatureEnum,VertexEnum); + UpdateInputsFromSolutionx(fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,tg,analysis_type, sub_analysis_type); + + if (!converged){ + if(verbose)_printf_("%s%i\n"," #unstable constraints = ",num_unstable_constraints); + if (num_unstable_constraints <= min_thermal_constraints)converged=1; + } + count++; + + if(converged==1)break; + } + + /*Free ressources: */ + MatFree(&Kgg_nopenalty); + VecFree(&pg_nopenalty); + VecFree(&tf); + VecFree(&tf_old); + delete solver_string; + + /*Assign output pointers: */ + *ptg=tg; + *pmelting_offset=melting_offset; +} Index: /issm/trunk/src/c/solvers/solvers.h =================================================================== --- /issm/trunk/src/c/solvers/solvers.h (revision 4029) +++ /issm/trunk/src/c/solvers/solvers.h (revision 4029) @@ -0,0 +1,19 @@ +/* + * solvers.h: + */ + +#ifndef SOLVERS_H_ +#define SOLVERS_H_ + +#include "../objects/objects.h" +#include "../io/io.h" + +struct OptArgs; +class FemModel; + +void solver_thermal_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem,int analysis_type,int sub_analysis_type); +void solver_diagnostic_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* fem,bool conserve_loads,int analysis_type,int sub_analysis_type); +void solver_linear(Vec* pug, FemModel* femmodel,int analysis_type,int sub_analysis_type); + + +#endif