Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 24079) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 24080) @@ -307,4 +307,5 @@ ./solutionsequences/solutionsequence_hydro_nonlinear.cpp\ ./solutionsequences/solutionsequence_shakti_nonlinear.cpp\ + ./solutionsequences/solutionsequence_glads_nonlinear.cpp\ ./cores/stressbalance_core.cpp\ ./solutionsequences/solutionsequence_stokescoupling_nonlinear.cpp\ Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 24079) +++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 24080) @@ -5,5 +5,5 @@ #include "../modules/modules.h" -#define AEPS 2.22e-15 +#define AEPS 2.2204460492503131E-015 /*Model processing*/ @@ -423,4 +423,18 @@ /*GlaDS specifics*/ +void HydrologyGlaDSAnalysis::SetChannelCrossSectionOld(FemModel* femmodel){/*{{{*/ + + bool ischannels; + femmodel->parameters->FindParam(&ischannels,HydrologyIschannelsEnum); + if(!ischannels) return; + + for(int i=0;iloads->Size();i++){ + if(femmodel->loads->GetEnum(i)==ChannelEnum){ + Channel* channel=(Channel*)femmodel->loads->GetObjectByOffset(i); + channel->SetChannelCrossSectionOld(); + } + } + +}/*}}}*/ void HydrologyGlaDSAnalysis::UpdateSheetThickness(FemModel* femmodel){/*{{{*/ @@ -458,5 +472,5 @@ Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input); Input* b_input = element->GetInput(BedEnum); _assert_(b_input); - Input* h_input = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input); + Input* hold_input = element->GetInput(HydrologySheetThicknessOldEnum);_assert_(hold_input); Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); @@ -472,5 +486,5 @@ vx_input->GetInputValue(&vx,gauss); vy_input->GetInputValue(&vy,gauss); - h_input->GetInputValue(&h_old,gauss); + hold_input->GetInputValue(&h_old,gauss); B_input->GetInputValue(&B,gauss); n_input->GetInputValue(&n,gauss); Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.h (revision 24079) +++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.h (revision 24080) @@ -32,4 +32,5 @@ /*Specific to GlaDS*/ + void SetChannelCrossSectionOld(FemModel* femmodel); void UpdateSheetThickness(FemModel* femmodel); void UpdateSheetThickness(Element* element); Index: /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp (revision 24079) +++ /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp (revision 24080) @@ -25,5 +25,5 @@ #define ALPHA_S 5./4. #define BETA_S 3./2. -#define AEPS 2.22e-15 +#define AEPS 2.2204460492503131E-015 /*Channel constructors and destructor*/ @@ -670,10 +670,5 @@ IssmDouble qc = - ks * pow(h,ALPHA_S) * pow(Ngrad,BETA_S-2.) * dphids; - /*Compute f factor*/ - IssmDouble fFactor = 0.; - if(this->S>0. || qc*dPw>0.){ - fFactor = lc * qc; - } - + /*Ice rate factor*/ A=pow(B,-n); @@ -682,16 +677,33 @@ IssmDouble N = phi_0 - phi; - IssmDouble alpha = 1./(rho_ice*L)*( - fabs(Qprime*pow(this->S,ALPHA_C-1.)*dphids) - + C*Qprime*pow(this->S,ALPHA_C-1.)*dPw - ) - 2./pow(n,n)*A*pow(fabs(N),n-1.)*N; - - IssmDouble beta = 1./(rho_ice*L)*( fabs(lc*qc*dphids) + C*fFactor*dPw ); - - /*Solve ODE*/ - this->S = ODE1(alpha,beta,this->Sold,dt,1); - - /*Make sure Area > 0*/ - if(this->S<0.) this->S = 0.; + bool converged = false; + int count = 0; + + while(!converged){ + + IssmDouble Snew = this->S; + + /*Compute f factor*/ + IssmDouble fFactor = 0.; + if(this->S>0. || qc*dPw>0.){ + fFactor = lc * qc; + } + + IssmDouble alpha = 1./(rho_ice*L)*( + fabs(Qprime*pow(Snew,ALPHA_C-1.)*dphids) + + C*Qprime*pow(Snew,ALPHA_C-1.)*dPw + ) - 2./pow(n,n)*A*pow(fabs(N),n-1.)*N; + + IssmDouble beta = 1./(rho_ice*L)*( fabs(lc*qc*dphids) + C*fFactor*dPw ); + + /*Solve ODE*/ + this->S = ODE1(alpha,beta,this->Sold,dt,2); + if(this->S<0.) this->S = 0.; + _assert_(!xIsNan(this->S)); + + count++; + + if(fabs((this->S - Snew)/(Snew+AEPS))<1e-8 || count>=10) converged = true; + } /*Clean up and return*/ Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 24079) +++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 24080) @@ -195,14 +195,15 @@ HydrologyGlaDSAnalysis* analysis = new HydrologyGlaDSAnalysis(); femmodel->SetCurrentConfiguration(HydrologyGlaDSAnalysisEnum); + + /*Set fields as old*/ InputDuplicatex(femmodel,HydraulicPotentialEnum,HydraulicPotentialOldEnum); - analysis->UpdateEffectivePressure(femmodel); - solutionsequence_shakti_nonlinear(femmodel); + InputDuplicatex(femmodel,HydrologySheetThicknessEnum,HydrologySheetThicknessOldEnum); + analysis->SetChannelCrossSectionOld(femmodel); + + /*Solve for new potential*/ + solutionsequence_glads_nonlinear(femmodel); if(VerboseSolution()) _printf0_(" updating effective pressure\n"); analysis->UpdateEffectivePressure(femmodel); - if(VerboseSolution()) _printf0_(" updating sheet thickness\n"); /*Uses N, so needs to come after*/ - analysis->UpdateSheetThickness(femmodel); - if(VerboseSolution()) _printf0_(" updating channels cross section\n"); - analysis->UpdateChannelCrossSection(femmodel); delete analysis; } Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24079) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24080) @@ -570,4 +570,5 @@ syn keyword cConstant HydrologySheetConductivityEnum syn keyword cConstant HydrologySheetThicknessEnum +syn keyword cConstant HydrologySheetThicknessOldEnum syn keyword cConstant HydrologyWaterVxEnum syn keyword cConstant HydrologyWaterVyEnum @@ -941,4 +942,5 @@ syn keyword cConstant ChannelEnum syn keyword cConstant ChannelAreaEnum +syn keyword cConstant ChannelAreaOldEnum syn keyword cConstant ClosedEnum syn keyword cConstant ColinearEnum @@ -1277,4 +1279,5 @@ syn keyword cType Cfsurfacesquare syn keyword cType Channel +syn keyword cType classes syn keyword cType Constraint syn keyword cType Constraints @@ -1283,6 +1286,6 @@ syn keyword cType ControlInput syn keyword cType Covertree +syn keyword cType DatasetInput syn keyword cType DataSetParam -syn keyword cType DatasetInput syn keyword cType Definition syn keyword cType DependentObject @@ -1297,6 +1300,6 @@ syn keyword cType ElementHook syn keyword cType ElementMatrix +syn keyword cType Elements syn keyword cType ElementVector -syn keyword cType Elements syn keyword cType ExponentialVariogram syn keyword cType ExternalResult @@ -1305,9 +1308,10 @@ syn keyword cType Friction syn keyword cType Gauss +syn keyword cType GaussianVariogram +syn keyword cType gaussobjects syn keyword cType GaussPenta syn keyword cType GaussSeg syn keyword cType GaussTetra syn keyword cType GaussTria -syn keyword cType GaussianVariogram syn keyword cType GenericExternalResult syn keyword cType GenericOption @@ -1317,4 +1321,5 @@ syn keyword cType Input syn keyword cType Inputs +syn keyword cType IntArrayInput syn keyword cType IntInput syn keyword cType IntMatParam @@ -1324,4 +1329,5 @@ syn keyword cType IssmDirectApplicInterface syn keyword cType IssmParallelDirectApplicInterface +syn keyword cType krigingobjects syn keyword cType Load syn keyword cType Loads @@ -1334,4 +1340,5 @@ syn keyword cType Matice syn keyword cType Matlitho +syn keyword cType matrixobjects syn keyword cType MatrixParam syn keyword cType Misfit @@ -1346,6 +1353,6 @@ syn keyword cType Observations syn keyword cType Option +syn keyword cType Options syn keyword cType OptionUtilities -syn keyword cType Options syn keyword cType Param syn keyword cType Parameters @@ -1360,10 +1367,10 @@ syn keyword cType Regionaloutput syn keyword cType Results +syn keyword cType Riftfront syn keyword cType RiftStruct -syn keyword cType Riftfront syn keyword cType Seg syn keyword cType SegInput +syn keyword cType Segment syn keyword cType SegRef -syn keyword cType Segment syn keyword cType SpcDynamic syn keyword cType SpcStatic @@ -1385,8 +1392,4 @@ syn keyword cType Vertex syn keyword cType Vertices -syn keyword cType classes -syn keyword cType gaussobjects -syn keyword cType krigingobjects -syn keyword cType matrixobjects syn keyword cType AdjointBalancethickness2Analysis syn keyword cType AdjointBalancethicknessAnalysis @@ -1407,6 +1410,6 @@ syn keyword cType FreeSurfaceBaseAnalysis syn keyword cType FreeSurfaceTopAnalysis +syn keyword cType GiaIvinsAnalysis syn keyword cType GLheightadvectionAnalysis -syn keyword cType GiaIvinsAnalysis syn keyword cType HydrologyDCEfficientAnalysis syn keyword cType HydrologyDCInefficientAnalysis Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24079) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24080) @@ -566,4 +566,5 @@ HydrologySheetConductivityEnum, HydrologySheetThicknessEnum, + HydrologySheetThicknessOldEnum, HydrologyWaterVxEnum, HydrologyWaterVyEnum, @@ -939,4 +940,5 @@ ChannelEnum, ChannelAreaEnum, + ChannelAreaOldEnum, ClosedEnum, ColinearEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24079) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24080) @@ -572,4 +572,5 @@ case HydrologySheetConductivityEnum : return "HydrologySheetConductivity"; case HydrologySheetThicknessEnum : return "HydrologySheetThickness"; + case HydrologySheetThicknessOldEnum : return "HydrologySheetThicknessOld"; case HydrologyWaterVxEnum : return "HydrologyWaterVx"; case HydrologyWaterVyEnum : return "HydrologyWaterVy"; @@ -943,4 +944,5 @@ case ChannelEnum : return "Channel"; case ChannelAreaEnum : return "ChannelArea"; + case ChannelAreaOldEnum : return "ChannelAreaOld"; case ClosedEnum : return "Closed"; case ColinearEnum : return "Colinear"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24079) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24080) @@ -584,4 +584,5 @@ else if (strcmp(name,"HydrologySheetConductivity")==0) return HydrologySheetConductivityEnum; else if (strcmp(name,"HydrologySheetThickness")==0) return HydrologySheetThicknessEnum; + else if (strcmp(name,"HydrologySheetThicknessOld")==0) return HydrologySheetThicknessOldEnum; else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum; else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"MeshVertexonbase")==0) return MeshVertexonbaseEnum; else if (strcmp(name,"MeshVertexonboundary")==0) return MeshVertexonboundaryEnum; - else if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"Misfit")==0) return MisfitEnum; + if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum; + else if (strcmp(name,"Misfit")==0) return MisfitEnum; else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum; else if (strcmp(name,"NewDamage")==0) return NewDamageEnum; @@ -751,9 +752,9 @@ else if (strcmp(name,"SmbTz")==0) return SmbTzEnum; else if (strcmp(name,"SmbV")==0) return SmbVEnum; - else if (strcmp(name,"SmbVmean")==0) return SmbVmeanEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"SmbVz")==0) return SmbVzEnum; + if (strcmp(name,"SmbVmean")==0) return SmbVmeanEnum; + else if (strcmp(name,"SmbVz")==0) return SmbVzEnum; else if (strcmp(name,"SmbW")==0) return SmbWEnum; else if (strcmp(name,"SmbWini")==0) return SmbWiniEnum; @@ -874,9 +875,9 @@ else if (strcmp(name,"Outputdefinition56")==0) return Outputdefinition56Enum; else if (strcmp(name,"Outputdefinition57")==0) return Outputdefinition57Enum; - else if (strcmp(name,"Outputdefinition58")==0) return Outputdefinition58Enum; else stage=8; } if(stage==8){ - if (strcmp(name,"Outputdefinition59")==0) return Outputdefinition59Enum; + if (strcmp(name,"Outputdefinition58")==0) return Outputdefinition58Enum; + else if (strcmp(name,"Outputdefinition59")==0) return Outputdefinition59Enum; else if (strcmp(name,"Outputdefinition5")==0) return Outputdefinition5Enum; else if (strcmp(name,"Outputdefinition60")==0) return Outputdefinition60Enum; @@ -964,4 +965,5 @@ else if (strcmp(name,"Channel")==0) return ChannelEnum; else if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum; + else if (strcmp(name,"ChannelAreaOld")==0) return ChannelAreaOldEnum; else if (strcmp(name,"Closed")==0) return ClosedEnum; else if (strcmp(name,"Colinear")==0) return ColinearEnum; @@ -996,10 +998,10 @@ else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum; else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum; - else if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum; - else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum; + if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum; + else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum; + else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum; else if (strcmp(name,"Element")==0) return ElementEnum; else if (strcmp(name,"ElementHook")==0) return ElementHookEnum; @@ -1119,10 +1121,10 @@ else if (strcmp(name,"Matice")==0) return MaticeEnum; else if (strcmp(name,"Matlitho")==0) return MatlithoEnum; - else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum; - else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum; + if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum; + else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum; + else if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum; else if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum; else if (strcmp(name,"MaxDivergence")==0) return MaxDivergenceEnum; @@ -1242,10 +1244,10 @@ else if (strcmp(name,"SubelementFriction2")==0) return SubelementFriction2Enum; else if (strcmp(name,"SubelementMelt1")==0) return SubelementMelt1Enum; - else if (strcmp(name,"SubelementMelt2")==0) return SubelementMelt2Enum; - else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum; + if (strcmp(name,"SubelementMelt2")==0) return SubelementMelt2Enum; + else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum; + else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum; else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum; else if (strcmp(name,"Tetra")==0) return TetraEnum; Index: /issm/trunk-jpl/src/c/solutionsequences/convergence.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/convergence.cpp (revision 24079) +++ /issm/trunk-jpl/src/c/solutionsequences/convergence.cpp (revision 24080) @@ -7,5 +7,5 @@ #include "../shared/shared.h" -void convergence(bool* pconverged, Matrix* Kff,Vector* pf,Vector* uf,Vector* old_uf,IssmDouble eps_res,IssmDouble eps_rel,IssmDouble eps_abs){/*{{{*/ +void convergence(bool* pconverged, Matrix* Kff,Vector* pf,Vector* uf,Vector* old_uf,IssmDouble eps_res,IssmDouble eps_rel,IssmDouble eps_abs){ /*output*/ @@ -139,5 +139,3 @@ /*assign output*/ *pconverged=converged; -}/*}}}*/ - - +} Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_glads_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_glads_nonlinear.cpp (revision 24080) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_glads_nonlinear.cpp (revision 24080) @@ -0,0 +1,96 @@ +/*!\file: solutionsequence_nonlinear.cpp + * \brief: core of a non-linear solution, using fixed-point method + */ + +#include "./solutionsequences.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" + +void solutionsequence_glads_nonlinear(FemModel* femmodel){ + + /*intermediary: */ + Matrix* Kff = NULL; + Matrix* Kfs = NULL; + Vector* ug = NULL; + Vector* uf = NULL; + Vector* old_uf = NULL; + Vector* pf = NULL; + Vector* df = NULL; + Vector* ys = NULL; + + /*parameters:*/ + int max_nonlinear_iterations; + IssmDouble eps_res,eps_rel,eps_abs; + HydrologyGlaDSAnalysis* analysis = new HydrologyGlaDSAnalysis(); + + /*Recover parameters (FIXME: from Stress balance for now :( )*/ + femmodel->parameters->FindParam(&max_nonlinear_iterations,StressbalanceMaxiterEnum); + femmodel->parameters->FindParam(&eps_res,StressbalanceRestolEnum); + femmodel->parameters->FindParam(&eps_rel,StressbalanceReltolEnum); + femmodel->parameters->FindParam(&eps_abs,StressbalanceAbstolEnum); + femmodel->UpdateConstraintsx(); + + int count_out=0; + bool converged_out=false; + int count_in=0; + bool converged_in=false; + + /*Start non-linear iteration using input velocity: */ + GetSolutionFromInputsx(&ug,femmodel); + Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters); + + while(!converged_out){ + + count_in=0; + converged_in=false; + + while(!converged_in){ + /*save pointer to old solution*/ + delete old_uf;old_uf=uf; + delete ug; + + SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); + CreateNodalConstraintsx(&ys,femmodel->nodes); + Reduceloadx(pf, Kfs, ys); delete Kfs; + femmodel->profiler->Start(SOLVER); + Solverx(&uf, Kff, pf, old_uf, df, femmodel->parameters); + femmodel->profiler->Stop(SOLVER); + Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete ys; + + convergence(&converged_in,Kff,pf,uf,old_uf,eps_res,eps_rel,eps_abs); delete Kff; delete pf; delete df; + InputUpdateFromSolutionx(femmodel,ug); + + /*Increase count: */ + count_in++; + if(count_in>=max_nonlinear_iterations && !converged_in){ + _printf0_(" maximum number of nonlinear iterations of inner loop (" << max_nonlinear_iterations << ") exceeded\n"); + converged_in = true; + } + } + if(VerboseConvergence()) _printf0_(setw(50) << left << " Inner loop converged in "<UpdateSheetThickness(femmodel); + analysis->UpdateChannelCrossSection(femmodel); + + /*Converged if inner loop converged in one solution*/ + //if(count_in==1) converged_out = true; + + /*Increase count: */ + count_out++; + if(count_out>=max_nonlinear_iterations && !converged_out){ + _printf0_(" maximum number of nonlinear iterations of outer loop (" << max_nonlinear_iterations << ") exceeded\n"); + converged_out = true; + } + } + + if(VerboseConvergence()) _printf0_("\n total number of iterations: " << count_out<<"x"<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,converged,ConvergedEnum); - for(;;){ /*save pointer to old solution*/ @@ -56,5 +53,4 @@ convergence(&converged,Kff,pf,uf,old_uf,eps_res,eps_rel,eps_abs); delete Kff; delete pf; delete df; - InputUpdateFromConstantx(femmodel,converged,ConvergedEnum); InputUpdateFromSolutionx(femmodel,ug); Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 24079) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 24080) @@ -15,4 +15,5 @@ void solutionsequence_hydro_nonlinear(FemModel* femmodel); void solutionsequence_shakti_nonlinear(FemModel* femmodel); +void solutionsequence_glads_nonlinear(FemModel* femmodel); void solutionsequence_nonlinear(FemModel* femmodel,bool conserve_loads); void solutionsequence_newton(FemModel* femmodel);