source: issm/trunk/src/m/solvers/solver_stokescoupling_nonlinear.m@ 6013

Last change on this file since 6013 was 6013, checked in by Mathieu Morlighem, 14 years ago

Added parameters to ReduceLoad

File size: 3.7 KB
Line 
1function femmodel=solver_couplingstokes_nonlinear(femmodel,conserve_loads)
2%SOLVER_COUPLINGSTOKES_NONLINEAR - core solver of coupling run
3%
4% Usage:
5% [femmodel]=solver_couplingstokes_nonlinear(femmodel,conserve_loads)
6
7 %Get parameters
8 kffpartitioning=femmodel.parameters.Kff;
9
10 %initialize solution vector
11 converged=0; count=1;
12
13 %First get ug=ug_horiz+ug_vert
14 femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
15 ug_horiz=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, femmodel.loads, femmodel.materials, femmodel.parameters);
16
17 femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum);
18 ug_vert=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, femmodel.loads, femmodel.materials, femmodel.parameters);
19
20 while(~converged),
21
22 %First compute the horizontal velocity
23 femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
24
25 %Start non-linear iteration using input velocity:
26 ug_horiz=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, femmodel.loads, femmodel.materials, femmodel.parameters);
27 uf_horiz=Reducevectorgtof( ug_horiz, femmodel.nodesets,femmodel.parameters);
28
29 %Update the solution to make sure that vx and vxold are similar
30 [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug_horiz);
31
32 %save pointer to old velocity
33 old_ug_horiz=ug_horiz;
34 old_uf_horiz=uf_horiz;
35
36 if kffpartitioning,
37 [K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
38 p_f = Reduceload( p_f, K_fs, femmodel.ys,femmodel.parameters);
39 else
40 [K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
41 [K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets,femmodel.parameters);
42 p_f = Reduceloadfromgtof( p_g, K_fs, femmodel.ys, femmodel.nodesets,femmodel.parameters);
43 end
44
45 uf_horiz=Solver(K_ff,p_f,old_uf_horiz,femmodel.parameters);
46 ug_horiz= Mergesolutionfromftog( uf_horiz, femmodel.ys, femmodel.nodesets,femmodel.parameters);
47
48 [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug_horiz);
49
50 %Figure out if convergence have been reached
51 converged=convergence(K_ff,p_f,uf_horiz,old_uf_horiz,femmodel.parameters);
52 converged=1;
53
54% %Then compute vertical velocity
55% femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum);
56% if kffpartitioning,
57% [K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
58% p_f = Reduceload( p_f, K_fs, femmodel.ys,femmodel.parameters);
59% else
60% [K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
61% [K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets,femmodel.parameters);
62% p_f = Reduceloadfromgtof( p_g, K_fs, femmodel.ys, femmodel.nodesets,femmodel.parameters);
63% end
64%
65% displaystring(femmodel.parameters.Verbose>1,'%s%g',' condition number of stiffness matrix: ',condest(K_ff));
66% uf_vert=Solver(K_ff,p_f,[],femmodel.parameters);
67%
68% ug_vert= Mergesolutionfromftog( uf_vert, femmodel.ys, femmodel.nodesets,femmodel.parameters);
69%
70% [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug_vert);
71
72 end
73end
Note: See TracBrowser for help on using the repository browser.