source: issm/oecreview/Archive/11319-11340/ISSM-11335-11336.diff

proj/ice/larour/issm-uci-clean/trunk-jpl/src/m/solvers/solver_newton.m

@@ +1 @@
+function femmodel=solver_newton(femmodel)
+%SOLVER_NEWTON - core solver of diagnostic run
+%
+%   Usage:
+%      [femmodel]=solver_newton(femmodel)
+
+        %Branch on partitioning schema requested
+        maxiter=femmodel.parameters.DiagnosticMaxiter;
+        configuration_type=femmodel.parameters.ConfigurationType;
+        [femmodel.nodes]=UpdateConstraints(femmodel.nodes,femmodel.constraints,femmodel.parameters);
+
+        %initialize solution vector
+        converged=0; count=1;
+
+        %Start non-linear iteration using input velocity: 
+        ug=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, femmodel.loads, femmodel.materials, femmodel.parameters);
+        uf=Reducevectorgtof( ug, femmodel.nodes,femmodel.parameters);
+
+        %Update the solution to make sure that vx and vxold are similar
+        [femmodel.elements femmodel.loads]=InputUpdateFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,double(converged),ConvergedEnum);
+        [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug);
+
+        while(~converged),
+
+                %save pointer to old velocity
+                old_ug=ug;
+                old_uf=uf;
+
+                %Solver forward model
+                [K_ff,K_fs,p_f,df,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+                ys=CreateNodalConstraints(femmodel.nodes,configuration_type);
+                p_f=Reduceload( p_f, K_fs, ys);
+                issmprintf(VerboseSolver(),'%s%g','      condition number of stiffness matrix: ',condest(K_ff));
+                uf=Solver(K_ff,p_f,old_uf,df,femmodel.parameters);
+                ug=Mergesolutionfromftog( uf, ys, femmodel.nodes,femmodel.parameters); 
+                [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug);
+
+                %Check convergence
+                converged=convergence(K_ff,p_f,uf,old_uf,femmodel.parameters);
+                if(converged==1) break; end
+                if(count>maxiter),
+                        issmprintf(true,'%s%i%s','      maximum number of iterations ',maxiter,' exceeded');
+                        break;
+                end
+
+                %Prepare next iteration using Newton's method
+                [K_ff,K_fs,p_f,df,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+                ys=CreateNodalConstraints(femmodel.nodes,configuration_type);
+                p_f=Reduceload( p_f, K_fs, ys);
+                pJf=p_f-K_ff*uf;
+                Jff=CreateJacobianMatrix(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,kmax);
+                duf=Solver(Jff,pJf,[],[],femmodel.parameters);
+                ys=0*ys;
+                uf=uf+duf;
+                ug=Mergesolutionfromftog(uf,ys,femmodel.nodes,femmodel.parameters); 
+                [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug);
+
+                %increase count
+                count=count+1;
+
+        end
+end