Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 14991) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 14992) @@ -331,7 +331,7 @@ ./analyses/EnvironmentFinalize.cpp\ ./analyses/ad_core.cpp\ - ./solvers/solver_linear.cpp\ - ./solvers/solver_nonlinear.cpp\ - ./solvers/solver_newton.cpp\ + ./solutionsequences/solutionsequence_linear.cpp\ + ./solutionsequences/solutionsequence_nonlinear.cpp\ + ./solutionsequences/solutionsequence_newton.cpp\ ./classes/objects/Options/Option.h\ ./classes/objects/Options/GenericOption.h\ @@ -390,5 +390,5 @@ ./analyses/thermal_core.cpp\ ./analyses/enthalpy_core.cpp\ - ./solvers/solver_thermal_nonlinear.cpp + ./solutionsequences/solutionsequence_thermal_nonlinear.cpp #}}} #Control sources {{{ @@ -444,5 +444,5 @@ ./analyses/adjointbalancethickness_core.cpp\ ./analyses/AdjointCorePointerFromSolutionEnum.cpp\ - ./solvers/solver_adjoint_linear.cpp + ./solutionsequences/solutionsequence_adjoint_linear.cpp #}}} @@ -464,5 +464,5 @@ ./modules/ModelProcessorx/HydrologyDCEfficient/CreateParametersHydrologyDCEfficient.cpp \ ./analyses/hydrology_core.cpp\ - ./solvers/solver_hydro_nonlinear.cpp + ./solutionsequences/solutionsequence_hydro_nonlinear.cpp #}}} #Diagnostic sources {{{ @@ -480,5 +480,5 @@ ./modules/ModelProcessorx/DiagnosticHutter/CreateLoadsDiagnosticHutter.cpp \ ./analyses/diagnostic_core.cpp\ - ./solvers/solver_stokescoupling_nonlinear.cpp + ./solutionsequences/solutionsequence_stokescoupling_nonlinear.cpp #}}} #Balanced sources {{{ Index: /issm/trunk-jpl/src/c/analyses/AdjointCorePointerFromSolutionEnum.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointCorePointerFromSolutionEnum.cpp (revision 14991) +++ /issm/trunk-jpl/src/c/analyses/AdjointCorePointerFromSolutionEnum.cpp (revision 14992) @@ -14,5 +14,5 @@ #include "../shared/Enum/Enum.h" #include "../modules/modules.h" -#include "../solvers/solvers.h" +#include "../solutionsequences/solutionsequences.h" void AdjointCorePointerFromSolutionEnum(void (**padjointcore)(FemModel*),int solutiontype){ Index: /issm/trunk-jpl/src/c/analyses/AnalysisConfiguration.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AnalysisConfiguration.cpp (revision 14991) +++ /issm/trunk-jpl/src/c/analyses/AnalysisConfiguration.cpp (revision 14992) @@ -14,5 +14,5 @@ #include "../shared/Enum/Enum.h" #include "../modules/modules.h" -#include "../solvers/solvers.h" +#include "../solutionsequences/solutionsequences.h" void AnalysisConfiguration(int** panalyses,int* pnumanalyses, int solutiontype){ Index: /issm/trunk-jpl/src/c/analyses/CorePointerFromSolutionEnum.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/CorePointerFromSolutionEnum.cpp (revision 14991) +++ /issm/trunk-jpl/src/c/analyses/CorePointerFromSolutionEnum.cpp (revision 14992) @@ -14,5 +14,5 @@ #include "../shared/Enum/Enum.h" #include "../modules/modules.h" -#include "../solvers/solvers.h" +#include "../solutionsequences/solutionsequences.h" void CorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype){ Index: /issm/trunk-jpl/src/c/analyses/WrapperCorePointerFromSolutionEnum.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/WrapperCorePointerFromSolutionEnum.cpp (revision 14991) +++ /issm/trunk-jpl/src/c/analyses/WrapperCorePointerFromSolutionEnum.cpp (revision 14992) @@ -13,5 +13,5 @@ #include "../shared/Enum/Enum.h" #include "../modules/modules.h" -#include "../solvers/solvers.h" +#include "../solutionsequences/solutionsequences.h" void WrapperCorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype,bool nodakotacore){ Index: /issm/trunk-jpl/src/c/analyses/ad_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ad_core.cpp (revision 14991) +++ /issm/trunk-jpl/src/c/analyses/ad_core.cpp (revision 14992) @@ -18,5 +18,5 @@ #include "../shared/Enum/Enum.h" #include "../modules/modules.h" -#include "../solvers/solvers.h" +#include "../solutionsequences/solutionsequences.h" /*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/adjointbalancethickness_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/adjointbalancethickness_core.cpp (revision 14991) +++ /issm/trunk-jpl/src/c/analyses/adjointbalancethickness_core.cpp (revision 14992) @@ -9,5 +9,5 @@ #include "../shared/Enum/Enum.h" #include "../modules/modules.h" -#include "../solvers/solvers.h" +#include "../solutionsequences/solutionsequences.h" void adjointbalancethickness_core(FemModel* femmodel){ @@ -22,5 +22,5 @@ if(VerboseSolution()) _pprintLine_(" computing thickness"); femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum); - solver_linear(femmodel); + solutionsequence_linear(femmodel); /*Call SurfaceAreax, because some it might be needed by PVector*/ @@ -30,5 +30,5 @@ if(VerboseSolution()) _pprintLine_(" computing adjoint"); femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum,AdjointBalancethicknessAnalysisEnum); - solver_adjoint_linear(femmodel); + solutionsequence_adjoint_linear(femmodel); /*Save results*/ Index: /issm/trunk-jpl/src/c/analyses/adjointdiagnostic_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/adjointdiagnostic_core.cpp (revision 14991) +++ /issm/trunk-jpl/src/c/analyses/adjointdiagnostic_core.cpp (revision 14992) @@ -9,5 +9,5 @@ #include "../shared/Enum/Enum.h" #include "../modules/modules.h" -#include "../solvers/solvers.h" +#include "../solutionsequences/solutionsequences.h" void adjointdiagnostic_core(FemModel* femmodel){ @@ -25,5 +25,5 @@ if(VerboseSolution()) _pprintLine_(" computing velocities"); femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum); - solver_nonlinear(femmodel,conserve_loads); + solutionsequence_nonlinear(femmodel,conserve_loads); /*Call SurfaceAreax, because some it might be needed by PVector*/ @@ -33,5 +33,5 @@ if(VerboseSolution()) _pprintLine_(" computing adjoint"); femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum,AdjointHorizAnalysisEnum); - solver_adjoint_linear(femmodel); + solutionsequence_adjoint_linear(femmodel); /*Save results*/ Index: /issm/trunk-jpl/src/c/analyses/balancethickness_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/balancethickness_core.cpp (revision 14991) +++ /issm/trunk-jpl/src/c/analyses/balancethickness_core.cpp (revision 14992) @@ -9,5 +9,5 @@ #include "../shared/Enum/Enum.h" #include "../modules/modules.h" -#include "../solvers/solvers.h" +#include "../solutionsequences/solutionsequences.h" void balancethickness_core(FemModel* femmodel){ @@ -23,5 +23,5 @@ if(VerboseSolution()) _pprintLine_("call computational core:"); - solver_linear(femmodel); + solutionsequence_linear(femmodel); if(save_results){ Index: /issm/trunk-jpl/src/c/analyses/bedslope_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/bedslope_core.cpp (revision 14991) +++ /issm/trunk-jpl/src/c/analyses/bedslope_core.cpp (revision 14992) @@ -9,5 +9,5 @@ #include "../shared/Enum/Enum.h" #include "../modules/modules.h" -#include "../solvers/solvers.h" +#include "../solutionsequences/solutionsequences.h" void bedslope_core(FemModel* femmodel){ @@ -23,7 +23,7 @@ /*Call on core computations: */ femmodel->SetCurrentConfiguration(BedSlopeAnalysisEnum,BedSlopeXAnalysisEnum); - solver_linear(femmodel); + solutionsequence_linear(femmodel); femmodel->SetCurrentConfiguration(BedSlopeAnalysisEnum,BedSlopeYAnalysisEnum); - solver_linear(femmodel); + solutionsequence_linear(femmodel); if(save_results){ Index: /issm/trunk-jpl/src/c/analyses/control_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/control_core.cpp (revision 14991) +++ /issm/trunk-jpl/src/c/analyses/control_core.cpp (revision 14992) @@ -9,5 +9,5 @@ #include "../shared/Enum/Enum.h" #include "../modules/modules.h" -#include "../solvers/solvers.h" +#include "../solutionsequences/solutionsequences.h" void control_core(FemModel* femmodel){ Index: /issm/trunk-jpl/src/c/analyses/controltao_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/controltao_core.cpp (revision 14991) +++ /issm/trunk-jpl/src/c/analyses/controltao_core.cpp (revision 14992) @@ -9,5 +9,5 @@ #include "../shared/Enum/Enum.h" #include "../modules/modules.h" -#include "../solvers/solvers.h" +#include "../solutionsequences/solutionsequences.h" #if defined (_HAVE_TAO_) && defined (_HAVE_PETSC_) && (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ > 1) Index: /issm/trunk-jpl/src/c/analyses/diagnostic_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/diagnostic_core.cpp (revision 14991) +++ /issm/trunk-jpl/src/c/analyses/diagnostic_core.cpp (revision 14992) @@ -9,5 +9,5 @@ #include "../shared/Enum/Enum.h" #include "../modules/modules.h" -#include "../solvers/solvers.h" +#include "../solutionsequences/solutionsequences.h" void diagnostic_core(FemModel* femmodel){ @@ -65,5 +65,5 @@ femmodel->SetCurrentConfiguration(DiagnosticHutterAnalysisEnum); - solver_linear(femmodel); + solutionsequence_linear(femmodel); if (ismacayealpattyn) ResetBoundaryConditions(femmodel,DiagnosticHorizAnalysisEnum); @@ -75,7 +75,7 @@ femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum); if(newton>0) - solver_newton(femmodel); + solutionsequence_newton(femmodel); else - solver_nonlinear(femmodel,conserve_loads); + solutionsequence_nonlinear(femmodel,conserve_loads); } @@ -83,5 +83,5 @@ if(VerboseSolution()) _pprintLine_(" computing coupling macayealpattyn and stokes velocities and pressure "); - solver_stokescoupling_nonlinear(femmodel,conserve_loads); + solutionsequence_stokescoupling_nonlinear(femmodel,conserve_loads); } @@ -90,5 +90,5 @@ if(VerboseSolution()) _pprintLine_(" computing vertical velocities"); femmodel->SetCurrentConfiguration(DiagnosticVertAnalysisEnum); - solver_linear(femmodel); + solutionsequence_linear(femmodel); } Index: /issm/trunk-jpl/src/c/analyses/enthalpy_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/enthalpy_core.cpp (revision 14991) +++ /issm/trunk-jpl/src/c/analyses/enthalpy_core.cpp (revision 14992) @@ -9,5 +9,5 @@ #include "../modules/modules.h" #include "../shared/io/io.h" -#include "../solvers/solvers.h" +#include "../solutionsequences/solutionsequences.h" void enthalpy_core(FemModel* femmodel){ @@ -21,5 +21,5 @@ if(VerboseSolution()) _pprintLine_(" computing enthalpy"); femmodel->SetCurrentConfiguration(EnthalpyAnalysisEnum); - solver_nonlinear(femmodel,true); + solutionsequence_nonlinear(femmodel,true); /*transfer enthalpy to enthalpy picard for the next step: */ Index: /issm/trunk-jpl/src/c/analyses/gia_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/gia_core.cpp (revision 14991) +++ /issm/trunk-jpl/src/c/analyses/gia_core.cpp (revision 14992) @@ -9,5 +9,5 @@ #include "../shared/Enum/Enum.h" #include "../modules/modules.h" -#include "../solvers/solvers.h" +#include "../solutionsequences/solutionsequences.h" void gia_core(FemModel* femmodel){ Index: /issm/trunk-jpl/src/c/analyses/gradient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/gradient_core.cpp (revision 14991) +++ /issm/trunk-jpl/src/c/analyses/gradient_core.cpp (revision 14992) @@ -9,5 +9,5 @@ #include "../shared/Enum/Enum.h" #include "../modules/modules.h" -#include "../solvers/solvers.h" +#include "../solutionsequences/solutionsequences.h" void gradient_core(FemModel* femmodel,int step,bool orthogonalize){ Index: /issm/trunk-jpl/src/c/analyses/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/hydrology_core.cpp (revision 14991) +++ /issm/trunk-jpl/src/c/analyses/hydrology_core.cpp (revision 14992) @@ -9,5 +9,5 @@ #include "../shared/Enum/Enum.h" #include "../modules/modules.h" -#include "../solvers/solvers.h" +#include "../solutionsequences/solutionsequences.h" void hydrology_core(FemModel* femmodel){ @@ -61,5 +61,5 @@ if(VerboseSolution()) _pprintLine_(" computing water column"); femmodel->SetCurrentConfiguration(HydrologyShreveAnalysisEnum); - solver_nonlinear(femmodel,modify_loads); + solutionsequence_nonlinear(femmodel,modify_loads); /*transfer water column thickness to old water column thickness: */ @@ -81,5 +81,5 @@ femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); if(VerboseSolution()) _pprintLine_(" computing water head"); - solver_hydro_nonlinear(femmodel); + solutionsequence_hydro_nonlinear(femmodel); if(save_results && ((i+1)%output_frequency==0 || (i+1)==nsteps)){ if(VerboseSolution()) _pprintLine_(" saving results "); Index: /issm/trunk-jpl/src/c/analyses/objectivefunction.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/objectivefunction.cpp (revision 14991) +++ /issm/trunk-jpl/src/c/analyses/objectivefunction.cpp (revision 14992) @@ -13,5 +13,5 @@ #include "../classes/objects/objects.h" #include "../shared/Enum/Enum.h" -#include "../solvers/solvers.h" +#include "../solutionsequences/solutionsequences.h" #include "../modules/modules.h" /*}}}*/ @@ -58,8 +58,8 @@ } else if (solution_type==DiagnosticSolutionEnum){ - solver_nonlinear(femmodel,conserve_loads); + solutionsequence_nonlinear(femmodel,conserve_loads); } else if (solution_type==BalancethicknessSolutionEnum){ - solver_linear(femmodel); + solutionsequence_linear(femmodel); } else if (solution_type==WeakBalancethicknessSolutionEnum){ Index: /issm/trunk-jpl/src/c/analyses/prognostic_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/prognostic_core.cpp (revision 14991) +++ /issm/trunk-jpl/src/c/analyses/prognostic_core.cpp (revision 14992) @@ -9,5 +9,5 @@ #include "../shared/Enum/Enum.h" #include "../modules/modules.h" -#include "../solvers/solvers.h" +#include "../solutionsequences/solutionsequences.h" void prognostic_core(FemModel* femmodel){ @@ -45,5 +45,5 @@ } if(VerboseSolution()) _pprintLine_(" call computational core"); - solver_linear(femmodel); + solutionsequence_linear(femmodel); if(save_results){ Index: /issm/trunk-jpl/src/c/analyses/steadystate_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/steadystate_core.cpp (revision 14991) +++ /issm/trunk-jpl/src/c/analyses/steadystate_core.cpp (revision 14992) @@ -14,5 +14,5 @@ #include "../shared/Enum/Enum.h" #include "../modules/modules.h" -#include "../solvers/solvers.h" +#include "../solutionsequences/solutionsequences.h" void steadystate_core(FemModel* femmodel){ Index: /issm/trunk-jpl/src/c/analyses/surfaceslope_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/surfaceslope_core.cpp (revision 14991) +++ /issm/trunk-jpl/src/c/analyses/surfaceslope_core.cpp (revision 14992) @@ -8,5 +8,5 @@ #include "../shared/io/io.h" #include "../shared/Enum/Enum.h" -#include "../solvers/solvers.h" +#include "../solutionsequences/solutionsequences.h" #include "../modules/modules.h" @@ -23,7 +23,7 @@ /*Call on core computations: */ femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SurfaceSlopeXAnalysisEnum); - solver_linear(femmodel); + solutionsequence_linear(femmodel); femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SurfaceSlopeYAnalysisEnum); - solver_linear(femmodel); + solutionsequence_linear(femmodel); if(save_results){ Index: /issm/trunk-jpl/src/c/analyses/thermal_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/thermal_core.cpp (revision 14991) +++ /issm/trunk-jpl/src/c/analyses/thermal_core.cpp (revision 14992) @@ -9,5 +9,5 @@ #include "../shared/Enum/Enum.h" #include "../modules/modules.h" -#include "../solvers/solvers.h" +#include "../solutionsequences/solutionsequences.h" void thermal_core(FemModel* femmodel){ @@ -28,9 +28,9 @@ if(VerboseSolution()) _pprintLine_(" computing temperatures"); femmodel->SetCurrentConfiguration(ThermalAnalysisEnum); - solver_thermal_nonlinear(femmodel); + solutionsequence_thermal_nonlinear(femmodel); if(VerboseSolution()) _pprintLine_(" computing melting"); femmodel->SetCurrentConfiguration(MeltingAnalysisEnum); - solver_linear(femmodel); + solutionsequence_linear(femmodel); if(save_results){ Index: /issm/trunk-jpl/src/c/analyses/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/transient_core.cpp (revision 14991) +++ /issm/trunk-jpl/src/c/analyses/transient_core.cpp (revision 14992) @@ -16,5 +16,5 @@ #include "../shared/Enum/Enum.h" #include "../modules/modules.h" -#include "../solvers/solvers.h" +#include "../solutionsequences/solutionsequences.h" void transient_core(FemModel* femmodel){ Index: /issm/trunk-jpl/src/c/main/issm.h =================================================================== --- /issm/trunk-jpl/src/c/main/issm.h (revision 14991) +++ /issm/trunk-jpl/src/c/main/issm.h (revision 14992) @@ -19,5 +19,5 @@ #include "../toolkits/toolkits.h" #include "../analyses/analyses.h" -#include "../solvers/solvers.h" +#include "../solutionsequences/solutionsequences.h" #include "../modules/modules.h" Index: /issm/trunk-jpl/src/c/solutionsequences/CMakeLists.txt =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/CMakeLists.txt (revision 14992) +++ /issm/trunk-jpl/src/c/solutionsequences/CMakeLists.txt (revision 14992) @@ -0,0 +1,16 @@ +# Subdirectories {{{ +# }}} +# Include Directory {{{ +include_directories(AFTER $ENV{ISSM_DIR}/src/c/solutionsequences) +# }}} +# CORE_SOURCES {{{ +set(CORE_SOURCES $ENV{ISSM_DIR}/src/c/solutionsequences/solutionsequence_linear.cpp + $ENV{ISSM_DIR}/src/c/solutionsequences/solutionsequence_newton.cpp + $ENV{ISSM_DIR}/src/c/solutionsequences/solutionsequence_nonlinear.cpp PARENT_SCOPE) +# }}} +# THERMAL_SOURCES {{{ +set(THERMAL_SOURCES $ENV{ISSM_DIR}/src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp PARENT_SCOPE) +# }}} +# DIAGNOSTIC_SOURCES {{{ +set(DIAGNOSTIC_SOURCES $ENV{ISSM_DIR}/src/c/solutionsequences/solutionsequence_stokescoupling_nonlinear.cpp PARENT_SCOPE) +# }}} Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_adjoint_linear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_adjoint_linear.cpp (revision 14992) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_adjoint_linear.cpp (revision 14992) @@ -0,0 +1,35 @@ +/*!\file: solutionsequence_adjoint_linear.cpp + * \brief: numerical core of linear solutions + */ + +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" + +void solutionsequence_adjoint_linear(FemModel* femmodel){ + /*This is axactly the same solutionsequence as solutionsequence_linear except that Reduceloadfromgtofx and Mergesolutionfromftogx + * use the flag "true" so that all spc are taken as 0*/ + + /*intermediary: */ + Matrix* Kff = NULL; + Matrix* Kfs = NULL; + Vector* ug = NULL; + Vector* uf = NULL; + Vector* pf = NULL; + Vector* df = NULL; + Vector* ys = NULL; + int configuration_type; + + /*Recover parameters: */ + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + femmodel->UpdateConstraintsx(); + + femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL); + CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); + Reduceloadx(pf, Kfs, ys,true); delete Kfs; //true means spc = 0 + Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); delete Kff; delete pf; delete df; + Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters,true); delete ys; //true means spc0 + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); + delete ug; delete uf; +} Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 14992) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 14992) @@ -0,0 +1,133 @@ +/*!\file: solutionsequence_hydro_nonlinear.cpp + * \brief: core of the hydro solution + */ + +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/io/io.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" + +void solutionsequence_hydro_nonlinear(FemModel* femmodel){ + + /*solution : */ + Vector* ug=NULL; + Vector* uf=NULL; + Vector* uf_old=NULL; + Vector* ys=NULL; + IssmDouble sediment_kmax; + + /*intermediary: */ + Matrix* Kff=NULL; + Matrix* Kfs=NULL; + Vector* pf=NULL; + Vector* df=NULL; + + bool converged,isefficientlayer; + int constraints_converged; + int num_unstable_constraints; + int count; + int hydro_maxiter; + + /*Recover parameters: */ + femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); + femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);//FIXME + femmodel->BasisIntegralsx(); + hydro_maxiter=100; + count=1; + converged=false; + + for(;;){ + + /*Computing the transfer term + + + */ + + + + /*First layer*/ + femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,true,ResetPenaltiesEnum); + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,false,ConvergedEnum); + femmodel->UpdateConstraintsx(); + femmodel->parameters->SetParam(HydrologySedimentEnum,HydrologyLayerEnum); + for(;;){ + femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df, &sediment_kmax); + CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCInefficientAnalysisEnum); + Reduceloadx(pf,Kfs,ys); delete Kfs; delete uf; + Solverx(&uf, Kff, pf,uf_old, df, femmodel->parameters); + delete uf_old; uf_old=uf->Duplicate(); + delete Kff; delete pf; delete ug; delete df; + Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters); delete ys; + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug); + + ConstraintsStatex(&constraints_converged,&num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + + if (!converged){ + if(VerboseConvergence()) _pprintLine_(" #unstable constraints = " << num_unstable_constraints); + if(num_unstable_constraints==0) converged = true; + if (count>=hydro_maxiter){ + _error_(" maximum number of iterations (" << hydro_maxiter << ") exceeded"); + } + } + count++; + + if(converged){ + femmodel->parameters->SetParam(sediment_kmax,HydrologySedimentKmaxEnum); + + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,converged,ConvergedEnum); + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug); + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,HydrologySedimentKmaxEnum); + break; + } + } + + /*Second layer*/ + if(!isefficientlayer) break; + femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum); + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,true,ResetPenaltiesEnum); + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,false,ConvergedEnum); + femmodel->UpdateConstraintsx(); + femmodel->parameters->SetParam(HydrologyEfficientEnum,HydrologyLayerEnum); + converged = false; + for(;;){ + femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df,NULL); + CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum); + Reduceloadx(pf,Kfs,ys); delete Kfs; delete uf; + Solverx(&uf, Kff, pf,uf_old, df, femmodel->parameters); + delete uf_old; uf_old=uf->Duplicate(); + delete Kff;delete pf; delete ug; delete df; + Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters); delete ys; + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug); + + ConstraintsStatex(&constraints_converged,&num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + + if (!converged){ + if(VerboseConvergence()) _pprintLine_(" #unstable constraints = " << num_unstable_constraints); + if(num_unstable_constraints==0) converged = true; + if (count>=hydro_maxiter){ + _error_(" maximum number of iterations (" << hydro_maxiter << ") exceeded"); + } + } + count++; + + if(converged){ + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,converged,ConvergedEnum); + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,MeltingOffsetEnum); + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug); + break; + } + } + + /*System convergence check*/ + break; + } + + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); + + /*Free ressources: */ + delete ug; + delete uf; + delete uf_old; +} Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_linear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_linear.cpp (revision 14992) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_linear.cpp (revision 14992) @@ -0,0 +1,39 @@ +/*!\file: solutionsequence_linear.cpp + * \brief: numerical core of linear solutions + */ + +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" + +void solutionsequence_linear(FemModel* femmodel){ + + /*intermediary: */ + Matrix* Kff = NULL; + Matrix* Kfs = NULL; + Vector* ug = NULL; + Vector* uf = NULL; + Vector* pf = NULL; + Vector* df = NULL; + Vector* ys = NULL; + int configuration_type; + + /*Recover parameters: */ + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + femmodel->UpdateConstraintsx(); + + femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL); + CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); + Reduceloadx(pf, Kfs, ys); delete Kfs; + Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); + delete Kff; delete pf; delete df; +//#ifdef _HAVE_ADOLC_ +// for (int i =0; isvector->M; ++i) { +// ADOLC_DUMP_MACRO(uf->svector->vector[i]); +// } +//#endif + Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys; + InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); + delete ug; +} Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_newton.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_newton.cpp (revision 14992) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_newton.cpp (revision 14992) @@ -0,0 +1,114 @@ +/*!\file: solutionsequence_nonlinear.cpp + * \brief: core of a non-linear solution, using fixed-point method + */ + +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/io/io.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" +#include "../analyses/analyses.h" + +void solutionsequence_newton(FemModel* femmodel){ + + /*intermediary: */ + bool converged; + int count,newton; + IssmDouble kmax; + Matrix* Kff = NULL; + Matrix* Kfs = NULL; + Matrix* Jff = NULL; + Vector* ug = NULL; + Vector* old_ug = NULL; + Vector* uf = NULL; + Vector* old_uf = NULL; + Vector* duf = NULL; + Vector* pf = NULL; + Vector* pJf = NULL; + Vector* df = NULL; + Vector* ys = NULL; + + /*parameters:*/ + int max_nonlinear_iterations; + int configuration_type; + + /*Recover parameters: */ + femmodel->parameters->FindParam(&max_nonlinear_iterations,DiagnosticMaxiterEnum); + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + femmodel->parameters->FindParam(&newton,DiagnosticIsnewtonEnum); + femmodel->UpdateConstraintsx(); + + count=1; + converged=false; + + /*Start non-linear iteration using input velocity: */ + GetSolutionFromInputsx(&ug,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + Reducevectorgtofx(&uf,ug,femmodel->nodes,femmodel->parameters); + + //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate) + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,converged,ConvergedEnum); + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug); + + for(;;){ + + delete old_ug;old_ug=ug; + delete old_uf;old_uf=uf; + + /*Solver forward model*/ + if(count==1 || newton==2){ + femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL); + CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); + Reduceloadx(pf,Kfs,ys);delete Kfs; + Solverx(&uf,Kff,pf,old_uf,df,femmodel->parameters);delete df; delete Kff; delete pf; + Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);delete ys; + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug); + delete old_ug;old_ug=ug; + delete old_uf;old_uf=uf; + } + uf=old_uf->Duplicate(); old_uf->Copy(uf); + + /*Prepare next iteration using Newton's method*/ + femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);delete df; + CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); + Reduceloadx(pf,Kfs,ys);delete Kfs; + + pJf=pf->Duplicate(); + Kff->MatMult(uf,pJf);// delete Kff); + pJf->Scale(-1.0); pJf->AXPY(pf,+1.0); //delete pf); + + femmodel->CreateJacobianMatrixx(&Jff,kmax); + Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters); delete Jff; delete pJf; + uf->AXPY(duf, 1.0); delete duf; + Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);delete ys; + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug); + + /*Check convergence*/ + convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); + delete Kff; delete pf; + if(converged==true){ + bool max_iteration_state=false; + int tempStep=1; + IssmDouble tempTime=1.0; + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1, MaxIterationConvergenceFlagEnum, max_iteration_state, tempStep, tempTime)); + break; + } + if(count>=max_nonlinear_iterations){ + _pprintLine_(" maximum number of Newton iterations (" << max_nonlinear_iterations << ") exceeded"); + bool max_iteration_state=true; + int tempStep=1; + IssmDouble tempTime=1.0; + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1, MaxIterationConvergenceFlagEnum, max_iteration_state, tempStep, tempTime)); + break; + } + + count++; + } + + if(VerboseConvergence()) _pprintLine_("\n total number of iterations: " << count-1); + + /*clean-up*/ + delete uf; + delete ug; + delete old_ug; + delete old_uf; +} Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp (revision 14992) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp (revision 14992) @@ -0,0 +1,118 @@ +/*!\file: solutionsequence_nonlinear.cpp + * \brief: core of a non-linear solution, using fixed-point method + */ + +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/io/io.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" +#include "../analyses/analyses.h" + +void solutionsequence_nonlinear(FemModel* femmodel,bool conserve_loads){ + + /*intermediary: */ + Matrix* Kff = NULL; + Matrix* Kfs = NULL; + Vector* ug = NULL; + Vector* old_ug = NULL; + Vector* uf = NULL; + Vector* old_uf = NULL; + Vector* pf = NULL; + Vector* df = NULL; + Vector* ys = NULL; + + Loads* savedloads=NULL; + bool converged; + int constraints_converged; + int num_unstable_constraints; + int count; + + /*parameters:*/ + int min_mechanical_constraints; + int max_nonlinear_iterations; + int configuration_type; + + /*Recover parameters: */ + femmodel->parameters->FindParam(&min_mechanical_constraints,DiagnosticRiftPenaltyThresholdEnum); + femmodel->parameters->FindParam(&max_nonlinear_iterations,DiagnosticMaxiterEnum); + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + femmodel->UpdateConstraintsx(); + + /*Were loads requested as output? : */ + if(conserve_loads){ + savedloads = static_cast(femmodel->loads->Copy()); + } + + count=1; + converged=false; + + /*Start non-linear iteration using input velocity: */ + GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); + Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters); + + //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate) + InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum); + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); + + for(;;){ + + //save pointer to old velocity + delete old_ug;old_ug=ug; + delete old_uf;old_uf=uf; + + femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL); + CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); + Reduceloadx(pf, Kfs, ys); delete Kfs; + Solverx(&uf, Kff, pf, old_uf, df, femmodel->parameters); + Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete ys; + + convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); delete Kff; delete pf; delete df; + InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum); + InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); + + ConstraintsStatex(&constraints_converged, &num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + if(VerboseConvergence()) _pprintLine_(" number of unstable constraints: " << num_unstable_constraints); + + //rift convergence + if (!constraints_converged) { + if (converged){ + if (num_unstable_constraints <= min_mechanical_constraints) converged=true; + else converged=false; + } + } + + /*Increase count: */ + count++; + if(converged==true){ + bool max_iteration_state=false; + int tempStep=1; + IssmDouble tempTime=1.0; + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1, MaxIterationConvergenceFlagEnum, max_iteration_state, tempStep, tempTime)); + break; + } + if(count>=max_nonlinear_iterations){ + _pprintLine_(" maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded"); + converged=true; + InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum); + InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); + bool max_iteration_state=true; + int tempStep=1; + IssmDouble tempTime=1.0; + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1, MaxIterationConvergenceFlagEnum, max_iteration_state, tempStep, tempTime)); + break; + } + } + + if(VerboseConvergence()) _pprintLine_("\n total number of iterations: " << count-1); + + /*clean-up*/ + if(conserve_loads){ + delete femmodel->loads; + femmodel->loads=savedloads; + } + delete uf; + delete ug; + delete old_ug; + delete old_uf; +} Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_stokescoupling_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_stokescoupling_nonlinear.cpp (revision 14992) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_stokescoupling_nonlinear.cpp (revision 14992) @@ -0,0 +1,101 @@ +/*!\file: solutionsequence_stokescoupling_nonlinear.cpp + * \brief: core of the coupling between stokes and macayealpattyn + */ + +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/Enum/Enum.h" +#include "../shared/io/io.h" +#include "../modules/modules.h" +#include "../analyses/analyses.h" + +void solutionsequence_stokescoupling_nonlinear(FemModel* femmodel,bool conserve_loads){ + + /*intermediary: */ + Matrix *Kff_horiz = NULL; + Matrix *Kfs_horiz = NULL; + Vector *ug_horiz = NULL; + Vector *uf_horiz = NULL; + Vector *old_uf_horiz = NULL; + Vector *pf_horiz = NULL; + Vector *df_horiz = NULL; + Matrix *Kff_vert = NULL; + Matrix *Kfs_vert = NULL; + Vector *ug_vert = NULL; + Vector *uf_vert = NULL; + Vector *pf_vert = NULL; + Vector *df_vert = NULL; + Vector *ys = NULL; + bool converged; + int count; + + /*parameters:*/ + int min_mechanical_constraints; + int max_nonlinear_iterations; + int configuration_type; + + /*Recover parameters: */ + femmodel->parameters->FindParam(&min_mechanical_constraints,DiagnosticRiftPenaltyThresholdEnum); + femmodel->parameters->FindParam(&max_nonlinear_iterations,DiagnosticMaxiterEnum); + femmodel->UpdateConstraintsx(); + + count=1; + converged=false; + + /*First get ug_horiz:*/ + femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum); + GetSolutionFromInputsx(&ug_horiz, femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters); + Reducevectorgtofx(&uf_horiz, ug_horiz, femmodel->nodes,femmodel->parameters); + + for(;;){ + + /*First diagnostic horiz:*/ + femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum); + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + + //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate) + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_horiz); + delete ug_horiz; + + //save pointer to old velocity + delete old_uf_horiz; old_uf_horiz=uf_horiz; + + /*solve: */ + femmodel->SystemMatricesx(&Kff_horiz, &Kfs_horiz, &pf_horiz, &df_horiz, NULL); + CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); + Reduceloadx(pf_horiz, Kfs_horiz, ys); delete Kfs_horiz; + Solverx(&uf_horiz, Kff_horiz, pf_horiz, old_uf_horiz, df_horiz,femmodel->parameters); + Mergesolutionfromftogx(&ug_horiz, uf_horiz,ys,femmodel->nodes,femmodel->parameters); delete ys; + InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_horiz); + + convergence(&converged,Kff_horiz,pf_horiz,uf_horiz,old_uf_horiz,femmodel->parameters); delete Kff_horiz; delete pf_horiz; delete df_horiz; + + /*Second compute vertical velocity: */ + femmodel->SetCurrentConfiguration(DiagnosticVertAnalysisEnum); + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + + /*solve: */ + femmodel->SystemMatricesx(&Kff_vert, &Kfs_vert, &pf_vert, &df_vert,NULL); + CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); + Reduceloadx(pf_vert, Kfs_vert, ys); delete Kfs_vert; + Solverx(&uf_vert, Kff_vert, pf_vert, NULL, df_vert,femmodel->parameters); delete Kff_vert; delete pf_vert; delete df_vert; + Mergesolutionfromftogx(&ug_vert, uf_vert,ys,femmodel->nodes,femmodel->parameters); + delete uf_vert; + delete ys; + InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_vert); + delete ug_vert; + + /*Increase count: */ + count++; + if(converged==true)break; + if(count>=max_nonlinear_iterations){ + _pprintLine_(" maximum number of iterations (" << max_nonlinear_iterations << ") exceeded"); + break; + } + } + + /*clean-up*/ + delete old_uf_horiz; + delete uf_horiz; + delete ug_horiz; +} Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp (revision 14992) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp (revision 14992) @@ -0,0 +1,83 @@ +/*!\file: solutionsequence_thermal_nonlinear.cpp + * \brief: core of the thermal solution + */ + +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/io/io.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" + +void solutionsequence_thermal_nonlinear(FemModel* femmodel){ + + /*solution : */ + Vector* tg=NULL; + Vector* tf=NULL; + Vector* tf_old=NULL; + Vector* ys=NULL; + IssmDouble melting_offset; + + /*intermediary: */ + Matrix* Kff=NULL; + Matrix* Kfs=NULL; + Vector* pf=NULL; + Vector* df=NULL; + + bool converged; + int constraints_converged; + int num_unstable_constraints; + int count; + int thermal_penalty_threshold; + int thermal_maxiter; + + /*parameters:*/ + int configuration_type; + + /*Recover parameters: */ + femmodel->parameters->FindParam(&thermal_penalty_threshold,ThermalPenaltyThresholdEnum); + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + femmodel->parameters->FindParam(&thermal_maxiter,ThermalMaxiterEnum); + + count=1; + converged=false; + + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,true,ResetPenaltiesEnum); + InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,false,ConvergedEnum); + femmodel->UpdateConstraintsx(); + + for(;;){ + + femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df, &melting_offset); + CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); + Reduceloadx(pf, Kfs, ys); delete Kfs; delete tf; + Solverx(&tf, Kff, pf,tf_old, df, femmodel->parameters); + delete tf_old; tf_old=tf->Duplicate(); + delete Kff;delete pf;delete tg; delete df; + Mergesolutionfromftogx(&tg, tf,ys,femmodel->nodes,femmodel->parameters); delete ys; + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,tg); + + ConstraintsStatex(&constraints_converged, &num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + + if (!converged){ + if(VerboseConvergence()) _pprintLine_(" #unstable constraints = " << num_unstable_constraints); + if (num_unstable_constraints <= thermal_penalty_threshold)converged=true; + if (count>=thermal_maxiter){ + converged=true; + _pprintLine_(" maximum number of iterations (" << thermal_maxiter << ") exceeded"); + } + } + count++; + + InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum); + + if(converged)break; + } + + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,tg); + InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,melting_offset,MeltingOffsetEnum); + + /*Free ressources: */ + delete tg; + delete tf; + delete tf_old; +} Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 14992) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 14992) @@ -0,0 +1,19 @@ +/* + * solutionsequences.h: + */ + +#ifndef _SOLUTION_SEQUENCES_H_ +#define _SOLUTION_SEQUENCES_H_ + +struct OptArgs; +class FemModel; + +void solutionsequence_thermal_nonlinear(FemModel* femmodel); +void solutionsequence_hydro_nonlinear(FemModel* femmodel); +void solutionsequence_nonlinear(FemModel* femmodel,bool conserve_loads); +void solutionsequence_newton(FemModel* femmodel); +void solutionsequence_stokescoupling_nonlinear(FemModel* femmodel,bool conserve_loads); +void solutionsequence_linear(FemModel* femmodel); +void solutionsequence_adjoint_linear(FemModel* femmodel); + +#endif