1 | function 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
|
---|
73 | end
|
---|