Changeset 11322

Timestamp:

02/06/12 07:45:50 (13 years ago)

Author:

Mathieu Morlighem

Message:

Added Newton's method for Pattyn, parallel only for now, to be improved

Location:

issm/trunk-jpl/src

Files:

• c/EnumDefinitions/EnumDefinitions.h (modified) (1 diff)
• c/Makefile.am (modified) (2 diffs)
• c/modules/CreateJacobianMatrixx (added)
• c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp (added)
• c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h (added)
• c/modules/EnumToStringx/EnumToStringx.cpp (modified) (1 diff)
• c/modules/ModelProcessorx/CreateParameters.cpp (modified) (1 diff)
• c/modules/StringToEnumx/StringToEnumx.cpp (modified) (1 diff)
• c/modules/modules.h (modified) (1 diff)
• c/objects/Elements/Element.h (modified) (1 diff)
• c/objects/Elements/Penta.cpp (modified) (2 diffs)
• c/objects/Elements/Penta.h (modified) (2 diffs)
• c/objects/Elements/Tria.cpp (modified) (4 diffs)
• c/objects/Elements/Tria.h (modified) (1 diff)
• c/objects/Materials/Matice.cpp (modified) (1 diff)
• c/objects/Materials/Matice.h (modified) (1 diff)
• c/objects/Numerics/ElementMatrix.cpp (modified) (2 diffs)
• c/objects/Numerics/ElementMatrix.h (modified) (1 diff)
• c/solutions/controltao_core.cpp (modified) (6 diffs)
• c/solutions/diagnostic_core.cpp (modified) (3 diffs)
• c/solvers/solver_newton.cpp (added)
• c/solvers/solver_nonlinear.cpp (modified) (1 diff)
• c/solvers/solvers.h (modified) (1 diff)
• m/classes/diagnostic.m (modified) (4 diffs)
• m/enum/DiagnosticIsnewtonEnum.m (added)
• m/enum/PentaP1ElementResultEnum.m (added)
• m/enum/PentaP1InputEnum.m (added)
• m/enum/PentaVertexElementResultEnum.m (deleted)
• m/enum/PentaVertexInputEnum.m (deleted)
• m/enum/TriaP1ElementResultEnum.m (added)
• m/enum/TriaP1InputEnum.m (added)
• m/enum/TriaVertexElementResultEnum.m (deleted)
• m/enum/TriaVertexInputEnum.m (deleted)
• m/model/loadresultsfromdisk.m (modified) (2 diffs)

issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h

@@ -28 +28 @@
         DiagnosticAbstolEnum,
         DiagnosticIcefrontEnum,
+        DiagnosticIsnewtonEnum,
         DiagnosticMaxiterEnum,
         DiagnosticNumRequestedOutputsEnum,

issm/trunk-jpl/src/c/Makefile.am

@@ -323 +323 @@
                                         ./modules/SystemMatricesx/SystemMatricesx.cpp\
                                         ./modules/SystemMatricesx/SystemMatricesx.h\
+                                        ./modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp\
+                                        ./modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h\
                                         ./modules/ConstraintsStatex/ConstraintsStatex.cpp\
                                         ./modules/ConstraintsStatex/ConstraintsStatex.h\
@@ -352 +354 @@
                                         ./solutions/CorePointerFromSolutionEnum.cpp\
                                         ./solvers/solver_linear.cpp\
-                                        ./solvers/solver_nonlinear.cpp
+                                        ./solvers/solver_nonlinear.cpp\
+                                        ./solvers/solver_newton.cpp
 #}}}
 #DAKOTA sources  {{{1

issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp

@@ -32 +32 @@
                 case DiagnosticAbstolEnum : return "DiagnosticAbstol";
                 case DiagnosticIcefrontEnum : return "DiagnosticIcefront";
+                case DiagnosticIsnewtonEnum : return "DiagnosticIsnewton";
                 case DiagnosticMaxiterEnum : return "DiagnosticMaxiter";
                 case DiagnosticNumRequestedOutputsEnum : return "DiagnosticNumRequestedOutputs";

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

@@ -40 +40 @@
         parameters->AddObject(iomodel->CopyConstantObject(DiagnosticReltolEnum));
         parameters->AddObject(iomodel->CopyConstantObject(DiagnosticAbstolEnum));
+        parameters->AddObject(iomodel->CopyConstantObject(DiagnosticIsnewtonEnum));
         parameters->AddObject(iomodel->CopyConstantObject(DiagnosticMaxiterEnum));
         parameters->AddObject(iomodel->CopyConstantObject(SteadystateReltolEnum));

issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp

@@ -30 +30 @@
         else if (strcmp(name,"DiagnosticAbstol")==0) return DiagnosticAbstolEnum;
         else if (strcmp(name,"DiagnosticIcefront")==0) return DiagnosticIcefrontEnum;
+        else if (strcmp(name,"DiagnosticIsnewton")==0) return DiagnosticIsnewtonEnum;
         else if (strcmp(name,"DiagnosticMaxiter")==0) return DiagnosticMaxiterEnum;
         else if (strcmp(name,"DiagnosticNumRequestedOutputs")==0) return DiagnosticNumRequestedOutputsEnum;

issm/trunk-jpl/src/c/modules/modules.h

@@ -108 +108 @@
 #include "./SurfaceAreax/SurfaceAreax.h"
 #include "./SystemMatricesx/SystemMatricesx.h"
+#include "./CreateJacobianMatrixx/CreateJacobianMatrixx.h"
 #include "./TimeAdaptx/TimeAdaptx.h"
 #include "./TriaSearchx/TriaSearchx.h"

issm/trunk-jpl/src/c/objects/Elements/Element.h

@@ -30 +30 @@
                 virtual void   CreateKMatrix(Mat Kff, Mat Kfs,Vec df)=0;
                 virtual void   CreatePVector(Vec pf)=0;
+                virtual void   CreateJacobianMatrix(Mat Jff)=0;
                 virtual void   GetSolutionFromInputs(Vec solution)=0;
                 virtual int    GetNodeIndex(Node* node)=0;

issm/trunk-jpl/src/c/objects/Elements/Penta.cpp

@@ -767 +767 @@
         /*clean up and return*/
         return pe;
+}
+/*}}}*/
+/*FUNCTION Penta::CreateJacobianMatrix{{{1*/
+void  Penta::CreateJacobianMatrix(Mat Jff){
+
+        /*retrieve parameters: */
+        ElementMatrix* Ke=NULL;
+        int analysis_type;
+        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+
+        /*Checks in debugging {{{2*/
+        _assert_(this->nodes && this->matice && this->matpar && this->verticalneighbors && this->parameters && this->inputs);
+        /*}}}*/
+
+        /*Skip if water element*/
+        if(IsOnWater()) return;
+
+        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
+        switch(analysis_type){
+#ifdef _HAVE_DIAGNOSTIC_
+                case DiagnosticHorizAnalysisEnum:
+                        Ke=CreateJacobianDiagnosticPattyn();
+                        break;
+#endif
+                default:
+                        _error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
+        }
+
+        /*Add to global matrix*/
+        if(Ke){
+                Ke->AddToGlobal(Jff);
+                delete Ke;
+        }
 }
 /*}}}*/
@@ -6224 +6257 @@
 }
 /*}}}*/
+/*FUNCTION Penta::CreateJacobianDiagnosticPattyn{{{1*/
+ElementMatrix* Penta::CreateJacobianDiagnosticPattyn(void){
+
+        /*Constants*/
+        const int    numdof=NDOF2*NUMVERTICES;
+
+        /*Intermediaries */
+        int        i,j,ig;
+        int        approximation;
+        double     xyz_list[NUMVERTICES][3];
+        double     Jdet;
+        double     eps1dotdphii,eps1dotdphij;
+        double     eps2dotdphii,eps2dotdphij;
+        double     mu_prime;
+        double     epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
+        double     eps1[3],eps2[3];
+        double     phi[NUMVERTICES];
+        double     dphi[3][NUMVERTICES];
+        GaussPenta *gauss=NULL;
+
+        /*Initialize Jacobian with regular Pattyn (first part of the Gateau derivative)*/
+        ElementMatrix* Ke=CreateKMatrixDiagnosticPattyn();
+
+        /*Retrieve all inputs and parameters*/
+        inputs->GetInputValue(&approximation,ApproximationEnum);
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+        Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
+        Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
+
+        /* Start  looping on the number of gaussian points: */
+        gauss=new GaussPenta(5,5);
+        for (ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+                GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
+
+                this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+                matice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
+                eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
+                eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];
+                eps1[2]=epsilon[3];                eps2[2]=epsilon[4];
+
+                for(i=0;i<6;i++){
+                        for(j=0;j<6;j++){
+                                eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i];
+                                eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j];
+                                eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i];
+                                eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j];
+
+                                Ke->values[12*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
+                                Ke->values[12*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
+                                Ke->values[12*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
+                                Ke->values[12*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
+                        }
+                }
+        }
+
+        /*Transform Coordinate System*/
+        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
+
+        /*Clean up and return*/
+        delete gauss;
+        return Ke;
+}
+/*}}}*/
 /*FUNCTION Penta::CreateKMatrixDiagnosticPattynFriction{{{1*/
 ElementMatrix* Penta::CreateKMatrixDiagnosticPattynFriction(void){

issm/trunk-jpl/src/c/objects/Elements/Penta.h

@@ -88 +88 @@
                 void   CreateKMatrix(Mat Kff, Mat Kfs,Vec df);
                 void   CreatePVector(Vec pf);
+                void   CreateJacobianMatrix(Mat Jff);
                 void   DeleteResults(void);
                 int    GetNodeIndex(Node* node);
@@ -234 +235 @@
                 ElementMatrix* CreateKMatrixDiagnosticVertVolume(void);
                 ElementMatrix* CreateKMatrixDiagnosticVertSurface(void);
+                ElementMatrix* CreateJacobianDiagnosticPattyn(void);
                 void           InputUpdateFromSolutionDiagnosticHoriz( double* solutiong);
                 void           InputUpdateFromSolutionDiagnosticMacAyeal( double* solutiong);

issm/trunk-jpl/src/c/objects/Elements/Tria.cpp

@@ -1883 +1883 @@
                                 name==InversionVyObsEnum ||
                                 name==FrictionCoefficientEnum ||
+                                name==MaterialsRheologyBbarEnum ||
                                 name==GradientEnum ||
                                 name==OldGradientEnum ||
@@ -4734 +4735 @@
 
         int doflist1[NUMVERTICES];
+        Input *input=NULL;
 
         /*Get out if this is not an element input*/
@@ -4742 +4744 @@
 
         /*Get input (either in element or material)*/
-        Input* input=inputs->GetInput(control_enum);
-        if(!input) _error_("Input %s not found in element",EnumToStringx(control_enum));
+        if(control_enum==MaterialsRheologyBbarEnum){
+                input=(Input*)matice->inputs->GetInput(control_enum); _assert_(input);
+        }
+        else{
+                input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
+        }
 
         /*Check that it is a ControlInput*/
@@ -4772 +4778 @@
         }
         new_input = new TriaP1Input(control_enum,values);
-
 
         if(control_enum==MaterialsRheologyBbarEnum){

issm/trunk-jpl/src/c/objects/Elements/Tria.h

@@ -84 +84 @@
                 void   CreateKMatrix(Mat Kff, Mat Kfs,Vec df);
                 void   CreatePVector(Vec pf);
+                void   CreateJacobianMatrix(Mat Jff){_error_("Not implemented yet");};
                 int    GetNodeIndex(Node* node);
                 int    Sid();

issm/trunk-jpl/src/c/objects/Materials/Matice.cpp

@@ -521 +521 @@
         /*Return: */
         *pviscosity_complement=viscosity_complement;
+}
+/*}}}*/
+/*FUNCTION Matice::GetViscosityDerivativeEpsSquare{{{1*/
+void  Matice::GetViscosityDerivativeEpsSquare(double* pmu_prime, double* epsilon){
+
+        /*output: */
+        double mu_prime;
+        double mu,n,eff2;
+
+        /*input strain rate: */
+        double exx,eyy,exy,exz,eyz;
+
+        /*Get visocisty and n*/
+        GetViscosity3d(&mu,epsilon);
+        n=GetN();
+
+        if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && 
+                                (epsilon[3]==0) && (epsilon[4]==0)){
+                mu_prime=0.5*pow((double)10,(double)14);
+        }
+        else{
+                /*Retrive strain rate components: */
+                exx=epsilon[0];
+                eyy=epsilon[1];
+                exy=epsilon[2];
+                exz=epsilon[3];
+                eyz=epsilon[4];
+                eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy + exz*exz + eyz*eyz;
+
+                mu_prime=(1-n)/(2*n) * mu/eff2;
+        }
+
+        /*Assign output pointers:*/
+        *pmu_prime=mu_prime;
 }
 /*}}}*/

issm/trunk-jpl/src/c/objects/Materials/Matice.h

@@ -68 +68 @@
                 void   GetViscosity3dStokes(double* pviscosity3d, double* epsilon);
                 void   GetViscosityComplement(double* pviscosity_complement, double* pepsilon);
+                void   GetViscosityDerivativeEpsSquare(double* pmu_prime, double* pepsilon);
                 double GetB();
                 double GetBbar();

issm/trunk-jpl/src/c/objects/Numerics/ElementMatrix.cpp

@@ -245 +245 @@
 
 /*ElementMatrix specific routines: */
-/*FUNCTION ElementMatrix::AddToGlobal{{{1*/
+/*FUNCTION ElementMatrix::AddToGlobal(Mat Kff, Mat Kfs){{{1*/
 void ElementMatrix::AddToGlobal(Mat Kff, Mat Kfs){
 
@@ -294 +294 @@
 }
 /*}}}*/
+/*FUNCTION ElementMatrix::AddToGlobal(Mat Jff){{{1*/
+void ElementMatrix::AddToGlobal(Mat Jff){
+
+        int i,j;
+        double* localvalues=NULL;
+
+        /*In debugging mode, check consistency (no NaN, and values not too big)*/
+        this->CheckConsistency();
+
+        if(this->dofsymmetrical){
+                /*only use row dofs to add values into global matrices: */
+
+                if(this->row_fsize){
+                        /*first, retrieve values that are in the f-set from the g-set values matrix: */
+                        localvalues=(double*)xmalloc(this->row_fsize*this->row_fsize*sizeof(double));
+                        for(i=0;i<this->row_fsize;i++){
+                                for(j=0;j<this->row_fsize;j++){
+                                        *(localvalues+this->row_fsize*i+j)=*(this->values+this->ncols*this->row_flocaldoflist[i]+this->row_flocaldoflist[j]);
+                                }
+                        }
+                        /*add local values into global  matrix, using the fglobaldoflist: */
+                        MatSetValues(Jff,this->row_fsize,this->row_fglobaldoflist,this->row_fsize,this->row_fglobaldoflist,(const double*)localvalues,ADD_VALUES);
+
+                        /*Free ressources:*/
+                        xfree((void**)&localvalues);
+                }
+
+        }
+        else{
+                _error_(" non dofsymmetrical matrix AddToGlobal routine not support yet!");
+        }
+
+}
+/*}}}*/
 /*FUNCTION ElementMatrix::CheckConsistency{{{1*/
 void ElementMatrix::CheckConsistency(void){

issm/trunk-jpl/src/c/objects/Numerics/ElementMatrix.h

@@ -59 +59 @@
                 /*ElementMatrix specific routines {{{1*/
                 void AddToGlobal(Mat Kff, Mat Kfs);
+                void AddToGlobal(Mat Jff);
                 void Echo(void);
                 void CheckConsistency(void);

issm/trunk-jpl/src/c/solutions/controltao_core.cpp

@@ -25 +25 @@
 
         /*TAO*/
-        int        ierr,numberofvertices;
+        int        ierr;
         int        num_controls;
         AppCtx     user;
@@ -45 +45 @@
 
         /*Initialize TAO*/
+        _printf_(VerboseControl(),"%s\n","   Initializing the Toolkit for Advanced Optimization (TAO)");
         TaoCreate(PETSC_COMM_WORLD,&tao);
         PetscOptionsSetValue("-tao_monitor","");
         TaoSetFromOptions(tao);
         TaoSetType(tao,"tao_blmvm");
+        //TaoSetType(tao,"tao_cg");
+        //TaoSetType(tao,"tao_lmvm");
 
         /*Prepare all TAO parameters*/
         TaoSetMaximumFunctionEvaluations(tao,50);
         TaoSetMaximumIterations(tao,10);
-        TaoSetTolerances(tao,10e-28,0.0000,0.0000,0.0000,10e-28);
-
-        numberofvertices = femmodel->vertices->NumberOfVertices();
-        XL=NewVec(numberofvertices); VecSet(XL,1.);
-        XU=NewVec(numberofvertices); VecSet(XU,200.);
+        TaoSetTolerances(tao,0.,0.,0.,0.,0.);
 
         GetVectorFromControlInputsx(&X, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
@@ -71 +70 @@
 
         /*Solver optimization problem*/
+        _printf_(VerboseControl(),"%s\n","   Starting optimization");
         TaoSolve(tao);
         TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);
@@ -96 +96 @@
 
         /*Set new variable*/
+        //VecView(X,PETSC_VIEWER_STDOUT_WORLD);
         SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
 
@@ -118 +119 @@
         /*Compute objective function*/
         CostFunctionx(fcn,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+        _printf_(true,"f(x)=%g\n",*fcn);
 
         /*Compute gradient*/
@@ -123 +125 @@
         VecCopy(gradient,G); VecFree(&gradient);
         VecScale(G,-1.);
+        //VecView(G,PETSC_VIEWER_STDOUT_WORLD);
 
         /*Clean-up and return*/

issm/trunk-jpl/src/c/solutions/diagnostic_core.cpp

@@ -16 +16 @@
 
         /*parameters: */
-        bool   dakota_analysis         = false;
-        int    dim                  = -1;
-        bool   ishutter             = false;
-        bool   ismacayealpattyn     = false;
-        bool   isstokes             = false;
-        bool   conserve_loads       = true;
-        bool   modify_loads         = true;
-        bool   control_analysis;
-        int    solution_type;
-        int    numoutputs         = 0;
-        int    *requested_outputs = NULL;
+        bool  dakota_analysis   = false;
+        int   dim               = -1;
+        bool  ishutter          = false;
+        bool  ismacayealpattyn  = false;
+        bool  isstokes          = false;
+        bool  isnewton          = false;
+        bool  conserve_loads    = true;
+        bool  modify_loads      = true;
+        bool  control_analysis;
+        int   solution_type;
+        int   numoutputs        = 0;
+        int  *requested_outputs = NULL;
 
         /* recover parameters:*/
@@ -33 +34 @@
         femmodel->parameters->FindParam(&ismacayealpattyn,FlowequationIsmacayealpattynEnum);
         femmodel->parameters->FindParam(&isstokes,FlowequationIsstokesEnum);
+        femmodel->parameters->FindParam(&isnewton,DiagnosticIsnewtonEnum);
         femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
         femmodel->parameters->FindParam(&control_analysis,InversionIscontrolEnum);
@@ -72 +74 @@
                 _printf_(VerboseSolution(),"%s\n","   computing velocities");
                 femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
-                solver_nonlinear(femmodel,modify_loads); 
+                if(isnewton)
+                 solver_newton(femmodel); 
+                else
+                 solver_nonlinear(femmodel,conserve_loads); 
         }
 

issm/trunk-jpl/src/c/solvers/solver_nonlinear.cpp

@@ -89 +89 @@
         }
 
+        _printf_(VerboseConvergence(),"\n   total number of iterations: %i\n",count-1);
+
         /*clean-up*/
         if(conserve_loads) delete loads;

issm/trunk-jpl/src/c/solvers/solvers.h

@@ -14 +14 @@
 void solver_thermal_nonlinear(FemModel* femmodel);
 void solver_nonlinear(FemModel* femmodel,bool conserve_loads);
+void solver_newton(FemModel* femmodel);
 void solver_stokescoupling_nonlinear(FemModel* femmodel,bool conserve_loads);
 void solver_linear(FemModel* femmodel);

issm/trunk-jpl/src/m/classes/diagnostic.m

@@ -12 +12 @@
                 reltol                   = 0;
                 abstol                   = 0;
+                isnewton                 = 0;
                 stokesreconditioning     = 0;
                 viscosity_overshoot      = 0;
@@ -82 +83 @@
                         checkfield(md,'diagnostic.reltol','size',[1 1]);
                         checkfield(md,'diagnostic.abstol','size',[1 1]);
+                        checkfield(md,'diagnostic.isnewton','numel',1,'values',[0 1]);
                         checkfield(md,'diagnostic.stokesreconditioning','size',[1 1],'NaN',1);
                         checkfield(md,'diagnostic.viscosity_overshoot','size',[1 1],'NaN',1);
@@ -127 +129 @@
                         fielddisplay(obj,'reltol','velocity relative convergence criterion, NaN -> not applied');
                         fielddisplay(obj,'abstol','velocity absolute convergence criterion, NaN -> not applied');
+                        fielddisplay(obj,'isnewton','Apply Newton''s method instead of a Picard fixed point method');
                         fielddisplay(obj,'maxiter','maximum number of nonlinear iterations');
                         fielddisplay(obj,'viscosity_overshoot','over-shooting constant new=new+C*(new-old)');
@@ -158 +161 @@
                         WriteData(fid,'object',obj,'fieldname','reltol','format','Double');
                         WriteData(fid,'object',obj,'fieldname','abstol','format','Double');
+                        WriteData(fid,'object',obj,'fieldname','isnewton','format','Boolean');
                         WriteData(fid,'object',obj,'fieldname','stokesreconditioning','format','Double');
                         WriteData(fid,'object',obj,'fieldname','viscosity_overshoot','format','Double');

issm/trunk-jpl/src/m/model/loadresultsfromdisk.m

@@ -10 +10 @@
         error('loadresultsfromdisk: error message.');
 end
+if ~exist(filename,'file'),
+        error(['binary file ' filename ' not found.']);
+end
 
 if ~md.qmu.isdakota,
@@ -20 +23 @@
         %load results onto model
         structure=parseresultsfromdisk(filename,~md.settings.io_gather);
+        if isempty(fieldnames(structure)),
+                error(['No result found in binary file ' filename '. Check for solution crash.']);
+        end
         md.results.(structure(1).SolutionType)=structure;