Changeset 5820

Timestamp:

09/15/10 09:52:40 (15 years ago)

Author:

seroussi

Message:

modifications in solver_stokescoupling

Location:

issm/trunk/src/m

Files:

• solutions/diagnostic_core.m (modified) (2 diffs)
• solvers/solver_stokescoupling_nonlinear.m (modified) (4 diffs)

issm/trunk/src/m/solutions/diagnostic_core.m

@@ -47 +47 @@
         end
                         
-        if ismacayealpattyn,
+        if xor(ismacayealpattyn,isstokes), %if macayealpattyn and stokes, we must use stokescoupling
 
                 displaystring(verbose,'\n%s',['computing horizontal velocities...']);
@@ -54 +54 @@
         end
         
+        if ismacayealpattyn & isstokes,
 
-        if dim==3 & ~isstokes,
+                displaystring(verbose,'\n%s',['computing coupling velocities for macayealpattyn and stokes...']);
+                femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
+                femmodel=solver_stokescoupling_nonlinear(femmodel,conserve_loads); 
+        end
+
+        if dim==3 & (ismacayealpattyn | ishutter),
         
                 displaystring(verbose,'\n%s',['computing vertical velocities...']);
                 femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum);
                 femmodel=solver_linear(femmodel);
-        end
-
-        if isstokes,
-
-                displaystring(verbose,'\n%s',['computing stokes velocities and pressure...']);
-                femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
-                femmodel=solver_diagnostic_nonlinear(femmodel,conserve_loads); 
         end
 

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

@@ -10 +10 @@
         %First get ug=ug_horiz+ug_vert
         femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
-        ug_horiz=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, loads, femmodel.materials, femmodel.parameters);
+        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, loads, femmodel.materials, femmodel.parameters);
-
-        ug=ug_horiz+ug_vert;
+        ug_vert=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, femmodel.loads, femmodel.materials, femmodel.parameters);
 
         while(~converged),
 
-                old_ug=ug;
-                
                 %First compute the horizontal velocity
                 femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
 
-                %keep a copy of loads for now
-                loads=femmodel.loads;
-
                 %Start non-linear iteration using input velocity: 
-                ug_horiz=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, loads, femmodel.materials, femmodel.parameters);
+                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
-                [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,ug_horiz);
+                [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug_horiz);
 
                 %save pointer to old velocity
@@ -38 +31 @@
                 old_uf_horiz=uf_horiz;
 
-                [K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters);
+                [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); 
@@ -47 +40 @@
                 ug_horiz= Mergesolutionfromftog( uf_horiz, femmodel.ys, femmodel.nodesets,femmodel.parameters); 
 
-                [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,ug_horiz);
+                [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug_horiz);
 
                 %Figure out if convergence have been reached
-                converged=convergence(K_ff,p_f,uf,old_uf,femmodel.parameters);
+                converged=convergence(K_ff,p_f,uf_horiz,old_uf_horiz,femmodel.parameters);
 
                 %Then compute vertical velocity
@@ -66 +59 @@
                 [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug_vert);
 
-                %Finally sum the two velocities and check the convergence
-                ug=ug_horiz+ug_vert;
-                converged=convergence([],[],ug,old_ug,femmodel.parameters);
-
         end