Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "PetscContactLineSearch.h"
11 : #include "FEProblem.h"
12 : #include "NonlinearSystem.h"
13 : #include "libmesh/petsc_nonlinear_solver.h"
14 : #include "libmesh/petsc_solver_exception.h"
15 : #include <petscdm.h>
16 :
17 : registerMooseObject("ContactApp", PetscContactLineSearch);
18 :
19 : InputParameters
20 86 : PetscContactLineSearch::validParams()
21 : {
22 86 : return ContactLineSearchBase::validParams();
23 : }
24 :
25 43 : PetscContactLineSearch::PetscContactLineSearch(const InputParameters & parameters)
26 43 : : ContactLineSearchBase(parameters)
27 : {
28 43 : _solver = dynamic_cast<PetscNonlinearSolver<Real> *>(
29 43 : _fe_problem.getNonlinearSystem(/*nl_sys_num=*/0).nonlinearSolver());
30 43 : if (!_solver)
31 0 : mooseError(
32 : "This line search operates only with Petsc, so Petsc must be your nonlinear solver.");
33 43 : }
34 :
35 : void
36 1182 : PetscContactLineSearch::lineSearch()
37 : {
38 1182 : PetscBool changed_y = PETSC_FALSE, changed_w = PETSC_FALSE;
39 : Vec X, F, Y, W, G, W1;
40 : SNESLineSearch line_search;
41 : PetscReal fnorm, xnorm, ynorm, gnorm;
42 : PetscReal ksp_rtol, ksp_abstol, ksp_dtol;
43 : PetscInt ksp_maxits;
44 : KSP ksp;
45 1182 : SNES snes = _solver->snes();
46 :
47 1182 : LibmeshPetscCall(SNESGetLineSearch(snes, &line_search));
48 1182 : LibmeshPetscCall(SNESLineSearchGetVecs(line_search, &X, &F, &Y, &W, &G));
49 1182 : LibmeshPetscCall(SNESLineSearchGetNorms(line_search, &xnorm, &fnorm, &ynorm));
50 1182 : LibmeshPetscCall(SNESLineSearchGetSNES(line_search, &snes));
51 1182 : LibmeshPetscCall(SNESLineSearchSetReason(line_search, SNES_LINESEARCH_SUCCEEDED));
52 1182 : LibmeshPetscCall(SNESGetKSP(snes, &ksp));
53 1182 : LibmeshPetscCall(KSPGetTolerances(ksp, &ksp_rtol, &ksp_abstol, &ksp_dtol, &ksp_maxits));
54 1182 : LibmeshPetscCall(VecDuplicate(W, &W1));
55 :
56 1182 : if (!_user_ksp_rtol_set)
57 : {
58 43 : _user_ksp_rtol = ksp_rtol;
59 43 : _user_ksp_rtol_set = true;
60 : }
61 :
62 1182 : ++_nl_its;
63 :
64 : /* precheck */
65 1182 : LibmeshPetscCall(SNESLineSearchPreCheck(line_search, X, Y, &changed_y));
66 :
67 : /* temporary update */
68 1182 : _contact_lambda = 1.;
69 1182 : LibmeshPetscCall(VecWAXPY(W, -_contact_lambda, Y, X));
70 :
71 : /* compute residual to determine whether contact state has changed since the last non-linear
72 : * residual evaluation */
73 : _current_contact_state.clear();
74 1182 : LibmeshPetscCall(SNESComputeFunction(snes, W, F));
75 : #if PETSC_VERSION_LESS_THAN(3, 25, 0)
76 : PetscBool domainerror;
77 : LibmeshPetscCall(SNESGetFunctionDomainError(snes, &domainerror));
78 : if (domainerror)
79 : LibmeshPetscCall(SNESLineSearchSetReason(line_search, SNES_LINESEARCH_FAILED_DOMAIN));
80 : #else
81 : SNESConvergedReason reason;
82 1182 : LibmeshPetscCall(SNESGetConvergedReason(snes, &reason));
83 1182 : if (reason == SNES_DIVERGED_FUNCTION_DOMAIN)
84 0 : LibmeshPetscCall(SNESLineSearchSetReason(line_search, SNES_LINESEARCH_FAILED_FUNCTION_DOMAIN));
85 : #endif
86 :
87 1182 : LibmeshPetscCall(VecNorm(F, NORM_2, &fnorm));
88 : std::set<dof_id_type> contact_state_stored = _current_contact_state;
89 : _current_contact_state.clear();
90 1182 : printContactInfo(contact_state_stored);
91 :
92 1182 : if (_affect_ltol)
93 : {
94 53 : if (contact_state_stored != _old_contact_state)
95 : {
96 19 : LibmeshPetscCall(KSPSetTolerances(ksp, _contact_ltol, ksp_abstol, ksp_dtol, ksp_maxits));
97 19 : _console << "Contact set changed since previous non-linear iteration!" << std::endl;
98 : }
99 : else
100 34 : LibmeshPetscCall(KSPSetTolerances(ksp, _user_ksp_rtol, ksp_abstol, ksp_dtol, ksp_maxits));
101 : }
102 :
103 : size_t ls_its = 0;
104 1739 : while (ls_its < _allowed_lambda_cuts)
105 : {
106 1470 : _contact_lambda *= 0.5;
107 : /* update */
108 1470 : LibmeshPetscCall(VecWAXPY(W1, -_contact_lambda, Y, X));
109 :
110 1470 : LibmeshPetscCall(SNESComputeFunction(snes, W1, G));
111 : #if PETSC_VERSION_LESS_THAN(3, 25, 0)
112 : LibmeshPetscCall(SNESGetFunctionDomainError(snes, &domainerror));
113 : if (domainerror)
114 : LibmeshPetscCall(SNESLineSearchSetReason(line_search, SNES_LINESEARCH_FAILED_DOMAIN));
115 : #else
116 1470 : LibmeshPetscCall(SNESGetConvergedReason(snes, &reason));
117 1470 : if (reason == SNES_DIVERGED_FUNCTION_DOMAIN)
118 0 : LibmeshPetscCall(
119 : SNESLineSearchSetReason(line_search, SNES_LINESEARCH_FAILED_FUNCTION_DOMAIN));
120 : #endif
121 :
122 1470 : LibmeshPetscCall(VecNorm(G, NORM_2, &gnorm));
123 1470 : if (gnorm < fnorm)
124 : {
125 557 : LibmeshPetscCall(VecCopy(G, F));
126 557 : LibmeshPetscCall(VecCopy(W1, W));
127 557 : fnorm = gnorm;
128 : contact_state_stored.swap(_current_contact_state);
129 : _current_contact_state.clear();
130 557 : printContactInfo(contact_state_stored);
131 557 : ++ls_its;
132 : }
133 : else
134 : break;
135 : }
136 :
137 1182 : LibmeshPetscCall(VecScale(Y, _contact_lambda));
138 : /* postcheck */
139 1182 : LibmeshPetscCall(SNESLineSearchPostCheck(line_search, X, Y, W, &changed_y, &changed_w));
140 :
141 1182 : if (changed_y)
142 0 : LibmeshPetscCall(VecWAXPY(W, -1., Y, X));
143 :
144 1182 : if (changed_w || changed_y)
145 : {
146 0 : LibmeshPetscCall(SNESComputeFunction(snes, W, F));
147 : #if PETSC_VERSION_LESS_THAN(3, 25, 0)
148 : LibmeshPetscCall(SNESGetFunctionDomainError(snes, &domainerror));
149 : if (domainerror)
150 : LibmeshPetscCall(SNESLineSearchSetReason(line_search, SNES_LINESEARCH_FAILED_DOMAIN));
151 : #else
152 0 : LibmeshPetscCall(SNESGetConvergedReason(snes, &reason));
153 0 : if (reason == SNES_DIVERGED_FUNCTION_DOMAIN)
154 0 : LibmeshPetscCall(
155 : SNESLineSearchSetReason(line_search, SNES_LINESEARCH_FAILED_FUNCTION_DOMAIN));
156 : #endif
157 : contact_state_stored.swap(_current_contact_state);
158 : _current_contact_state.clear();
159 0 : printContactInfo(contact_state_stored);
160 : }
161 :
162 : /* copy the solution over */
163 1182 : LibmeshPetscCall(VecCopy(W, X));
164 :
165 1182 : LibmeshPetscCall(SNESLineSearchComputeNorms(line_search));
166 :
167 1182 : LibmeshPetscCall(VecDestroy(&W1));
168 :
169 : _old_contact_state = std::move(contact_state_stored);
170 1182 : }
|
