Changeset 1415

Timestamp:

07/29/09 15:40:59 (15 years ago)

Author:

Mathieu Morlighem

Message:

improved mechanical residue

File:

• issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m (modified) (3 diffs)

issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m

@@ -49 +49 @@
                 [soln(count).u_f]=Solver(K_ff,p_f,[],m.parameters);
 
-                %Get residue
-                residu=norm(K_ff*soln(count).u_f-p_f,2)/norm(p_f,2);
-
                 %Merge back to g set
                 [soln(count).u_g]= Mergesolutionfromftog( soln(count).u_f, m.Gmn, m.ys, m.nodesets ); 
-
-                if (count>2),
-                        soln(count-1).u_f=NaN;
-                end
 
                 %Deal with penalty loads
@@ -65 +58 @@
                 [loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( m.elements,m.nodes, loads, m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
 
-        %   Figure out if convergence is reached.
+                %Figure out if convergence have been reached
+
+                %Residue criterion
+                solver_residue=norm(K_ff*soln(count).u_f-p_f,2)/norm(p_f,2);
+                displaystring(m.parameters.debug,'%-53s%g','      convergence criterion: norm(KU-F)/norm(F)',solver_residue);
+
+                %Force equilibrium
+                if count>2,
+                        mechanical_residue=norm(K_ff*soln(count-1).u_f-p_f,2)/norm(p_f,2);
+                else
+                        mechanical_residue=NaN;
+                end
+                displaystring(m.parameters.debug,'%-53s%g%s','      convergence criterion: norm(KUold-F)/norm(F)',mechanical_residue*100,' %');
+
+                %Relative criterion
                 dug=soln(count).u_g-soln(count-1).u_g; 
-                nduinf=norm(dug,inf)*m.parameters.yts; 
                 ndu=norm(dug,2); 
                 nu=norm(soln(count-1).u_g,2); 
-
-                %Residue criterion
-                displaystring(m.parameters.debug,'%s%g','      convergence criterion: norm(KU-F)/norm(F)=',residu);
-
-                %Relative criterion
                 if (ndu/nu<=m.parameters.eps_rel),
-                        displaystring(m.parameters.debug,'%s%g%s%g','      convergence criterion: norm(du)/norm(u)=',ndu/nu,' < ',m.parameters.eps_rel);
+                        displaystring(m.parameters.debug,'%-53s%g%s%g%s','      convergence criterion: norm(du)/norm(u)',ndu/nu*100,' < ',m.parameters.eps_rel*100,' %');
                         converged=1;
                 else
-                        displaystring(m.parameters.debug,'%s%g%s%g','      convergence criterion: norm(du)/norm(u)=',ndu/nu,' > ',m.parameters.eps_rel);
+                        displaystring(m.parameters.debug,'%-53s%g%s%g%s','      convergence criterion: norm(du)/norm(u)',ndu/nu*100,' > ',m.parameters.eps_rel*100,' %');
                         converged=0;
                 end
 
                 %Absolute criterion
+                nduinf=norm(dug,inf)*m.parameters.yts; 
                 if ~isnan(m.parameters.eps_abs),
                         if (nduinf<=m.parameters.eps_abs),
-                                displaystring(m.parameters.debug,'%s%g%s%g','      convergence criterion: max(du)=',nduinf,' < ',m.parameters.eps_abs);
+                                displaystring(m.parameters.debug,'%-53s%g%s%g%s','      convergence criterion: max(du)',nduinf,' < ',m.parameters.eps_abs,' m/yr');
                         else
-                                displaystring(m.parameters.debug,'%s%g%s%g','      convergence criterion: max(du)=',nduinf,' > ',m.parameters.eps_abs);
+                                displaystring(m.parameters.debug,'%-53s%g%s%g%s','      convergence criterion: max(du)',nduinf,' > ',m.parameters.eps_abs,' m/yr');
                                 converged=0;
                         end
@@ -97 +99 @@
                         converged=0;
                 end
-
-                %if (count>2),
-                %       soln(count-1).u_g=NaN;
-                %end
         end
                         
         %some cleanup
-        [soln(count).u_f]=NaN;
-
         u_g=soln(count).u_g;
+        clear soln
 
         %more output might be needed, when running in cielocontrol.m
         nout=max(nargout,1)-1;
         if nout==2,
-                inputs=add(inputs,'velocity',soln(count).u_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes);
+                inputs=add(inputs,'velocity',u_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes);
                 m.parameters.kflag=1; m.parameters.pflag=0; 
                 [K_gg, p_g]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);