source: issm/branches/trunk-jpl-damage/src/c/solvers/solver_stokescoupling_nonlinear.cpp@ 11427

Last change on this file since 11427 was 11284, checked in by seroussi, 13 years ago

minor ConvergedEnum should always be a bool
File size: 4.1 KB
Line 
1/*!\file: solver_stokescoupling_nonlinear.cpp
2 * \brief: core of the coupling between stokes and macayealpattyn
3 */
4
5#include "../toolkits/toolkits.h"
6#include "../objects/objects.h"
7#include "../EnumDefinitions/EnumDefinitions.h"
8#include "../io/io.h"
9#include "../modules/modules.h"
10#include "../solutions/solutions.h"
11#include "./solvers.h"
12
13void solver_stokescoupling_nonlinear(FemModel* femmodel,bool conserve_loads){
14
15 /*intermediary: */
16 Mat Kff_horiz = NULL, Kfs_horiz = NULL;
17 Vec ug_horiz = NULL, uf_horiz = NULL, old_uf_horiz = NULL;
18 Vec pf_horiz = NULL;
19 Vec df_horiz = NULL;
20 Mat Kff_vert = NULL, Kfs_vert = NULL;
21 Vec ug_vert = NULL, uf_vert = NULL;
22 Vec pf_vert = NULL;
23 Vec df_vert = NULL;
24 Vec ys = NULL;
25 bool converged;
26 int constraints_converged;
27 int num_unstable_constraints;
28 int count;
29
30 /*parameters:*/
31 int min_mechanical_constraints;
32 int max_nonlinear_iterations;
33 int configuration_type;
34
35 /*Recover parameters: */
36 femmodel->parameters->FindParam(&min_mechanical_constraints,DiagnosticRiftPenaltyThresholdEnum);
37 femmodel->parameters->FindParam(&max_nonlinear_iterations,DiagnosticMaxiterEnum);
38 UpdateConstraintsx(femmodel->nodes,femmodel->constraints,femmodel->parameters);
39
40 count=1;
41 converged=false;
42
43 /*First get ug_horiz:*/
44 femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
45 GetSolutionFromInputsx(&ug_horiz, femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters);
46 Reducevectorgtofx(&uf_horiz, ug_horiz, femmodel->nodes,femmodel->parameters);
47
48 for(;;){
49
50 /*First diagnostic horiz:*/
51 femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
52 femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
53
54 //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate)
55 InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_horiz);
56 VecFree(&ug_horiz);
57
58 //save pointer to old velocity
59 VecFree(&old_uf_horiz);old_uf_horiz=uf_horiz;
60
61 /*solve: */
62 SystemMatricesx(&Kff_horiz, &Kfs_horiz, &pf_horiz, &df_horiz, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
63 CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
64 Reduceloadx(pf_horiz, Kfs_horiz, ys); MatFree(&Kfs_horiz);
65 Solverx(&uf_horiz, Kff_horiz, pf_horiz, old_uf_horiz, df_horiz,femmodel->parameters);
66 Mergesolutionfromftogx(&ug_horiz, uf_horiz,ys,femmodel->nodes,femmodel->parameters); VecFree(&ys);
67 InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_horiz);
68
69 convergence(&converged,Kff_horiz,pf_horiz,uf_horiz,old_uf_horiz,femmodel->parameters); MatFree(&Kff_horiz);VecFree(&pf_horiz); VecFree(&df_horiz);
70
71 /*Second compute vertical velocity: */
72 femmodel->SetCurrentConfiguration(DiagnosticVertAnalysisEnum);
73 femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
74
75 /*solve: */
76 SystemMatricesx(&Kff_vert, &Kfs_vert, &pf_vert, &df_vert,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
77 CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
78 Reduceloadx(pf_vert, Kfs_vert, ys); MatFree(&Kfs_vert);
79 Solverx(&uf_vert, Kff_vert, pf_vert, NULL, df_vert,femmodel->parameters); MatFree(&Kff_vert); VecFree(&pf_vert); VecFree(&df_vert);
80 Mergesolutionfromftogx(&ug_vert, uf_vert,ys,femmodel->nodes,femmodel->parameters);VecFree(&uf_vert); VecFree(&ys);
81 InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_vert);
82 VecFree(&ug_vert); VecFree(&uf_vert);
83
84 /*Increase count: */
85 count++;
86 if(converged==true)break;
87 if(count>=max_nonlinear_iterations){
88 _printf_(true," maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations);
89 break;
90 }
91 }
92
93 /*clean-up*/
94 VecFree(&old_uf_horiz);
95 VecFree(&uf_horiz);
96 VecFree(&ug_horiz);
97 VecFree(&ys);
98}
Note: See TracBrowser for help on using the repository browser.