Changeset 6023

Timestamp:

09/24/10 11:19:25 (14 years ago)

Author:

seroussi

Message:

finished solver for stokes coupling in serial

File:

• issm/trunk/src/m/solvers/solver_stokescoupling_nonlinear.m (modified) (4 diffs)

issm/trunk/src/m/solvers/solver_stokescoupling_nonlinear.m

@@ -14 +14 @@
         femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
         ug_horiz=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, femmodel.loads, femmodel.materials, femmodel.parameters);
-
-        femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum);
-        ug_vert=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, femmodel.loads, femmodel.materials, femmodel.parameters);
+        uf_horiz=Reducevectorgtof( ug_horiz, femmodel.nodesets,femmodel.parameters);
 
         while(~converged),
@@ -22 +20 @@
                 %First compute the horizontal velocity
                 femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
-
-                %Start non-linear iteration using input velocity: 
-                ug_horiz=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, femmodel.loads, femmodel.materials, femmodel.parameters);
-                uf_horiz=Reducevectorgtof( ug_horiz, femmodel.nodesets,femmodel.parameters);
 
                 %Update the solution to make sure that vx and vxold are similar
@@ -35 +29 @@
 
                 if kffpartitioning,
-                        [K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
-                        p_f = Reduceload( p_f, K_fs, femmodel.ys,femmodel.parameters);
+                        [K_gg_horiz,K_ff_horiz,K_fs_horiz,p_g_horiz,p_f_horiz,kmax_horiz]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+                        p_f_horiz = Reduceload( p_f_horiz, K_fs_horiz, femmodel.ys,femmodel.parameters);
                 else
-                        [K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
-                        [K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets,femmodel.parameters); 
-                        p_f = Reduceloadfromgtof( p_g, K_fs, femmodel.ys, femmodel.nodesets,femmodel.parameters);
+                        [K_gg_horiz,K_ff_horiz,K_fs_horiz,p_g_horiz,p_f_horiz,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+                        [K_ff_horiz, K_fs_horiz] = Reducematrixfromgtof( K_gg_horiz, femmodel.nodesets,femmodel.parameters); 
+                        p_f_horiz = Reduceloadfromgtof( p_g_horiz, K_fs_horiz, femmodel.ys, femmodel.nodesets,femmodel.parameters);
                 end
 
-                uf_horiz=Solver(K_ff,p_f,old_uf_horiz,femmodel.parameters);
+                uf_horiz=Solver(K_ff_horiz,p_f_horiz,old_uf_horiz,femmodel.parameters);
                 ug_horiz= Mergesolutionfromftog( uf_horiz, femmodel.ys, femmodel.nodesets,femmodel.parameters); 
 
@@ -49 +43 @@
 
                 %Figure out if convergence have been reached
-                converged=convergence(K_ff,p_f,uf_horiz,old_uf_horiz,femmodel.parameters);
-                converged=1;
+                converged=convergence(K_ff_horiz,p_f_horiz,uf_horiz,old_uf_horiz,femmodel.parameters);
 
-%               %Then compute vertical velocity
-%               femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum);
-%       if kffpartitioning,
-%               [K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
-%               p_f = Reduceload( p_f, K_fs, femmodel.ys,femmodel.parameters);
-%       else
-%               [K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
-%               [K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets,femmodel.parameters); 
-%               p_f = Reduceloadfromgtof( p_g, K_fs, femmodel.ys, femmodel.nodesets,femmodel.parameters);
-%       end
-%
-%               displaystring(femmodel.parameters.Verbose>1,'%s%g','      condition number of stiffness matrix: ',condest(K_ff));
-%               uf_vert=Solver(K_ff,p_f,[],femmodel.parameters);
-%
-%               ug_vert= Mergesolutionfromftog( uf_vert, femmodel.ys, femmodel.nodesets,femmodel.parameters); 
-%
-%               [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug_vert);
+                %Then compute vertical velocity
+                femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum);
+
+                if kffpartitioning,
+                        [K_gg_vert,K_ff_vert,K_fs_vert,p_g_vert,p_f_vert,kmax_vert]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+                        p_f_vert = Reduceload( p_f_vert, K_fs_vert, femmodel.ys,femmodel.parameters);
+                else
+                        [K_gg_vert,K_ff_vert,K_fs_vert,p_g_vert,p_f_vert,kmax_vert]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+                        [K_ff_vert, K_fs_vert] = Reducematrixfromgtof( K_gg_vert, femmodel.nodesets,femmodel.parameters); 
+                        p_f_vert = Reduceloadfromgtof( p_g_vert, K_fs_vert, femmodel.ys, femmodel.nodesets,femmodel.parameters);
+                end
+
+                uf_vert=Solver(K_ff_vert,p_f_vert,[],femmodel.parameters);
+
+                ug_vert= Mergesolutionfromftog( uf_vert, femmodel.ys, femmodel.nodesets,femmodel.parameters); 
+
+                [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug_vert);
 
         end