1 /* Siconos is a program dedicated to modeling, simulation and control
2 * of non smooth dynamical systems.
3 *
4 * Copyright 2021 INRIA.
5 *
6 * Licensed under the Apache License, Version 2.0 (the "License");
7 * you may not use this file except in compliance with the License.
8 * You may obtain a copy of the License at
9 *
10 * http://www.apache.org/licenses/LICENSE-2.0
11 *
12 * Unless required by applicable law or agreed to in writing, software
13 * distributed under the License is distributed on an "AS IS" BASIS,
14 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15 * See the License for the specific language governing permissions and
16 * limitations under the License.
17 */
18 #include <assert.h> // for assert
19 #include <stdio.h> // for NULL, fprintf, stderr
20 #include <stdlib.h> // for exit, EXIT_FAILURE
21 #include "Friction_cst.h" // for SICONOS_GLOBAL_FRICTION_3D...
22 #include "GlobalFrictionContactProblem.h" // for GlobalFrictionContactProblem
23 #include "NonSmoothDrivers.h" // for gfc3d_driver
24 #include "NumericsFwd.h" // for SolverOptions, GlobalFrict...
25 #include "SolverOptions.h" // for SolverOptions, solver_opti...
26 #include "siconos_debug.h" // for DEBUG_EXPR
27 #include "gfc3d_Solvers.h" // for gfc3d_ACLMFixedPoint, gfc3...
28 #include "numerics_verbose.h" // for numerics_printf_verbose
29 #include "gfc3d_balancing.h"
30 #include "gfc3d_compute_error.h"
31 #include "SiconosBlas.h" // for cblas_dcopy, cblas_dscal
32
33 #ifdef DEBUG_MESSAGES
34 #include "NumericsVector.h"
35 #include "NumericsMatrix.h"
36 #endif
37
38 const char* const SICONOS_GLOBAL_FRICTION_3D_NSGS_WR_STR = "GFC3D_NSGS_WR";
39 const char* const SICONOS_GLOBAL_FRICTION_3D_NSN_AC_WR_STR = "GFC3D_NSN_AC_WR";
40 const char* const SICONOS_GLOBAL_FRICTION_3D_NSGSV_WR_STR = "GFC3D_NSGSV_WR";
41 const char* const SICONOS_GLOBAL_FRICTION_3D_PROX_WR_STR = "GFC3D_PROX_WR";
42 const char* const SICONOS_GLOBAL_FRICTION_3D_DSFP_WR_STR = "GFC3D_DSFP_WR";
43 const char* const SICONOS_GLOBAL_FRICTION_3D_TFP_WR_STR = "GFC3D_TFP_WR";
44 const char* const SICONOS_GLOBAL_FRICTION_3D_NSGS_STR = "GFC3D_NSGS";
45 const char* const SICONOS_GLOBAL_FRICTION_3D_NSN_AC_STR = "GFC3D_NSN_AC";
46 const char* const SICONOS_GLOBAL_FRICTION_3D_GAMS_PATH_STR = "GFC3D_GAMS_PATH";
47 const char* const SICONOS_GLOBAL_FRICTION_3D_GAMS_PATHVI_STR = "GFC3D_GAMS_PATHVI";
48 const char* const SICONOS_GLOBAL_FRICTION_3D_VI_EG_STR = "GFC3D_VI_EG";
49 const char* const SICONOS_GLOBAL_FRICTION_3D_ACLMFP_STR = "GFC3D_ACLMFP";
50 const char* const SICONOS_GLOBAL_FRICTION_3D_VI_FPP_STR = "GFC3D_VI_FPP";
51 const char* const SICONOS_GLOBAL_FRICTION_3D_ADMM_WR_STR = "GFC3D_ADMM_WR";
52
53
gfc3d_balancing_check_drift(GlobalFrictionContactProblem * balanced_problem,GlobalFrictionContactProblem * problem,double * reaction,double * velocity,double * globalVelocity,SolverOptions * options)54 static int gfc3d_balancing_check_drift(GlobalFrictionContactProblem* balanced_problem,
55 GlobalFrictionContactProblem* problem,
56 double *reaction, double *velocity,
57 double* globalVelocity, SolverOptions* options)
58 {
59 if(options->iparam[SICONOS_FRICTION_3D_IPARAM_RESCALING]>0)
60 {
61 size_t nc = problem->numberOfContacts;
62 size_t n = problem->M->size0;
63 size_t m = 3 * nc;
64
65 double norm_b = cblas_dnrm2(m, balanced_problem->b, 1);
66 double norm_q = cblas_dnrm2(n, balanced_problem->q, 1);
67 double error_balancing = 1e24;
68 double tolerance = options->dparam[SICONOS_DPARAM_TOL];
69 gfc3d_compute_error(balanced_problem, reaction, velocity, globalVelocity,
70 tolerance, options,
71 norm_q, norm_b, &error_balancing);
72
73 /* Come back to original variables */
74 gfc3d_balancing_back_to_original_variables(balanced_problem, options,
75 reaction, velocity, globalVelocity);
76
77 norm_b = cblas_dnrm2(m, problem->b, 1);
78 norm_q = cblas_dnrm2(n, problem->q, 1);
79 double error =0.0;
80 gfc3d_compute_error(problem, reaction, velocity, globalVelocity,
81 tolerance, options,
82 norm_q, norm_b, &error);
83
84 numerics_printf_verbose(0,"error with balancing = %8.4e", error_balancing);
85 numerics_printf_verbose(0,"error with original = %8.4e", error);
86 }
87 //else continue
88
89 return 0;
90
91 }
92
93
gfc3d_driver(GlobalFrictionContactProblem * problem,double * reaction,double * velocity,double * globalVelocity,SolverOptions * options)94 int gfc3d_driver(GlobalFrictionContactProblem* problem, double *reaction, double *velocity,
95 double* globalVelocity, SolverOptions* options)
96 {
97 assert(options->isSet);
98 DEBUG_EXPR(NV_display(globalVelocity,problem_ori->M->size0););
99 if(verbose > 0)
100 solver_options_print(options);
101
102 /* Solver name */
103 /* const char* const name = options->solverName;*/
104
105
106
107 int info = -1 ;
108
109 if(problem->dimension != 3)
110 numerics_error("gfc3d_driver", "Dimension of the problem : problem-> dimension is not compatible or is not set");
111
112 /* if there is no contact, we compute directly the global velocity as M^{-1}q */
113 int m = problem->H->size1;
114 if(m ==0)
115 {
116 numerics_printf_verbose(1,"---- GFC3D - DRIVER . No contact case. Direct computation of global velocity");
117 globalFrictionContact_computeGlobalVelocity(problem, reaction, globalVelocity);
118 return 0;
119 }
120
121 /* Non Smooth Gauss Seidel (NSGS) */
122 switch(options->solverId)
123 {
124 case SICONOS_GLOBAL_FRICTION_3D_NSGS_WR:
125 {
126
127 numerics_printf_verbose(1," ========================== Call NSGS_WR solver with reformulation into Friction-Contact 3D problem ==========================\n");
128 gfc3d_nsgs_wr(problem, reaction, velocity, globalVelocity, &info, options);
129 break;
130
131 }
132 case SICONOS_GLOBAL_FRICTION_3D_NSGSV_WR:
133 {
134
135 numerics_printf_verbose(1," ========================== Call NSGSV_WR solver with reformulation into Friction-Contact 3D problem ==========================\n");
136 gfc3d_nsgs_velocity_wr(problem, reaction, velocity, globalVelocity, &info, options);
137 break;
138 }
139 case SICONOS_GLOBAL_FRICTION_3D_NSN_AC_WR:
140 {
141
142 numerics_printf_verbose(1," ========================== Call NSN_AC_WR solver with reformulation into Friction-Contact 3D problem ==========================\n");
143 gfc3d_nonsmooth_Newton_AlartCurnier_wr(problem, reaction, velocity, globalVelocity, &info, options);
144 break;
145
146 }
147 case SICONOS_GLOBAL_FRICTION_3D_PROX_WR:
148 {
149
150 numerics_printf_verbose(1," ========================== Call PROX_WR solver with reformulation into Friction-Contact 3D problem ==========================\n");
151 gfc3d_proximal_wr(problem, reaction, velocity, globalVelocity, &info, options);
152 break;
153
154 }
155 case SICONOS_GLOBAL_FRICTION_3D_DSFP_WR:
156 {
157
158 numerics_printf_verbose(1," ========================== Call DSFP_WR solver with reformulation into Friction-Contact 3D problem ==========================\n");
159 gfc3d_DeSaxceFixedPoint_wr(problem, reaction, velocity, globalVelocity, &info, options);
160 break;
161
162 }
163 case SICONOS_GLOBAL_FRICTION_3D_TFP_WR:
164 {
165
166 numerics_printf_verbose(1," ========================== Call TFP_WR solver with reformulation into Friction-Contact 3D problem ==========================\n");
167 gfc3d_TrescaFixedPoint_wr(problem, reaction, velocity, globalVelocity, &info, options);
168 break;
169
170 }
171 case SICONOS_GLOBAL_FRICTION_3D_NSGS:
172 {
173 gfc3d_nsgs(problem, reaction, velocity, globalVelocity,
174 &info, options);
175 break;
176
177 }
178 case SICONOS_GLOBAL_FRICTION_3D_NSN_AC:
179 {
180 /* Balancing */
181 /* here, the balancing is done outside the solver */
182 /* therfore the solver does not account for the possible drift of error measure between
183 the balanced problem and the original one */
184
185 GlobalFrictionContactProblem* balanced_problem = gfc3d_balancing_problem(problem,options);
186 gfc3d_balancing_go_to_balanced_variables(balanced_problem, options,
187 reaction, velocity, globalVelocity);
188 /* Call the solver with balanced data */
189 gfc3d_nonsmooth_Newton_AlartCurnier(balanced_problem, reaction, velocity,
190 globalVelocity, &info, options);
191
192 /* check if the drift is large */
193 // int info_check_drift =
194 gfc3d_balancing_check_drift(balanced_problem,problem, reaction, velocity, globalVelocity,
195 options);
196
197 problem = gfc3d_balancing_free(problem, options);
198
199 break;
200
201 }
202 case SICONOS_GLOBAL_FRICTION_3D_GAMS_PATH:
203 {
204 numerics_printf_verbose(1," ========================== Call PATH solver via GAMS for an AVI Friction-Contact 3D problem ==========================\n");
205 gfc3d_AVI_gams_path(problem, reaction, velocity, &info, options);
206 break;
207 }
208 case SICONOS_GLOBAL_FRICTION_3D_GAMS_PATHVI:
209 {
210 numerics_printf_verbose(1," ========================== Call PATHVI solver via GAMS for an AVI Friction-Contact 3D problem ==========================\n");
211 gfc3d_AVI_gams_pathvi(problem, reaction, globalVelocity, &info, options);
212 break;
213 }
214 case SICONOS_GLOBAL_FRICTION_3D_VI_FPP:
215 {
216 gfc3d_VI_FixedPointProjection(problem, reaction, velocity,
217 globalVelocity, &info, options);
218 break;
219
220 }
221 case SICONOS_GLOBAL_FRICTION_3D_VI_EG:
222 {
223 gfc3d_VI_ExtraGradient(problem, reaction, velocity,
224 globalVelocity, &info, options);
225 break;
226
227 }
228 case SICONOS_GLOBAL_FRICTION_3D_ACLMFP:
229 {
230 gfc3d_ACLMFixedPoint(problem, reaction, velocity,
231 globalVelocity, &info, options);
232 break;
233
234 }
235 case SICONOS_GLOBAL_FRICTION_3D_ADMM:
236 {
237 /* globalFrictionContact_rescaling(problem, 1.0/1.512808e-04, 1.0/1.407230e+01, 1.0); */
238 gfc3d_ADMM(problem, reaction, velocity,
239 globalVelocity, &info, options);
240 break;
241
242 }
243 case SICONOS_GLOBAL_FRICTION_3D_ADMM_WR:
244 {
245
246 numerics_printf_verbose(1," ========================== Call NSGS_WR solver with reformulation into Friction-Contact 3D problem ==========================\n");
247 gfc3d_admm_wr(problem, reaction, velocity, globalVelocity, &info, options);
248 break;
249
250 }
251 case SICONOS_GLOBAL_FRICTION_3D_IPM:
252 {
253 gfc3d_IPM(problem, reaction, velocity,
254 globalVelocity, &info, options);
255 break;
256
257 }
258 default:
259 {
260 fprintf(stderr, "Numerics, gfc3d_driver failed. Unknown solver %d.\n", options->solverId);
261 exit(EXIT_FAILURE);
262
263 }
264 }
265
266 return info;
267
268 }
269
gfc3d_checkTrivialCaseGlobal(int n,double * q,double * velocity,double * reaction,double * globalVelocity,SolverOptions * options)270 int gfc3d_checkTrivialCaseGlobal(int n, double* q, double* velocity, double* reaction, double * globalVelocity, SolverOptions* options)
271 {
272 /* norm of vector q */
273 /* double qs = cblas_dnrm2( n , q , 1 ); */
274 /* int i; */
275 int info = -1;
276 /* if( qs <= DBL_EPSILON ) */
277 /* { */
278 /* // q norm equal to zero (less than DBL_EPSILON) */
279 /* // -> trivial solution: reaction = 0 and velocity = q */
280 /* for( i = 0 ; i < n ; ++i ) */
281 /* { */
282 /* velocity[i] = q[i]; */
283 /* reaction[i] = 0.; */
284 /* } */
285 /* iparam[2] = 0; */
286 /* iparam[4]= 0; */
287 /* dparam[1] = 0.0; */
288 /* dparam[3] = 0.0; */
289 /* info = 0; */
290 /* if(iparam[1]>0) */
291 /* printf("fc3d driver, norm(q) = 0, trivial solution reaction = 0, velocity = q.\n"); */
292 /* } */
293 return info;
294 }
295