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

Index: /proj/ice/larour/issm-uci-clean/trunk-jpl/src/m/solvers/solver_newton.m
===================================================================
--- /proj/ice/larour/issm-uci-clean/trunk-jpl/src/m/solvers/solver_newton.m     (revision 0)
+++ /proj/ice/larour/issm-uci-clean/trunk-jpl/src/m/solvers/solver_newton.m     (revision 11336)
@@ -0,0 +1,62 @@
+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