GCC Code Coverage Report

Directory: ./
File: PETScInterface/SNESInterface.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 63 79 79.7%
Branches: 39 117 33.3%
Line Branch Exec Source
1 // ______ ______ _ _ _____ ______
2 // | ____| ____| | (_)/ ____| | ____|
3 // | |__ | |__ | | _| (___ ___| |__
4 // | __| | __| | | | |\___ \ / __| __|
5 // | | | |____| |____| |____) | (__| |____
6 // |_| |______|______|_|_____/ \___|______|
7 // Finite Elements for Life Sciences and Engineering
8 //
9 // License: LGL2.1 License
10 // FELiScE default license: LICENSE in root folder
11 //
12 // Main authors: Vicente Mataix Ferrandiz
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Core/mpiInfo.hpp"
21 #include "Solver/linearProblem.hpp"
22 #include "PETScInterface/SNESInterface.hpp"
23
24 namespace felisce
25 {
26 namespace { // anonymous
27 /**
28 * @brief This function is used to monitor evolution of the convergence.
29 * Is is intended to be passed to SNESMonitorSet()
30 * @param[in] snes SNES object. Not explicitly used but required by Petsc prototype.
31 * @param[in] its Index of iteration.
32 * @param[in] norm Current L^2 norm.
33 * @param[in] contextAsVoid Optional user-defined function context. In our case is is a pointer to LinearProblem object.
34 */
35 20911 PetscErrorCode snesViewerLinearProblem(
36 SNES snes,
37 PetscInt its,
38 PetscReal norm,
39 void* contextAsVoid
40 )
41 {
42 PetscFunctionBegin;
43
44
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 20911 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 20911 FEL_ASSERT(contextAsVoid);
45
46 20911 const LinearProblem* p_linear_problem = static_cast<const LinearProblem*>(contextAsVoid);
47
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 20911 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 20911 FEL_ASSERT(p_linear_problem);
48
49 PetscReal evaluation_state_Min, evaluation_state_Max;
50
51 20911 const PetscVector& r_vector = p_linear_problem->evaluationState();
52
1/2
✓ Branch 1 taken 20911 times.
✗ Branch 2 not taken.
 20911 r_vector.min(&evaluation_state_Min);
53
1/2
✓ Branch 1 taken 20911 times.
✗ Branch 2 not taken.
 20911 r_vector.max(&evaluation_state_Max);
54
55 20911 PetscObject snes_object = reinterpret_cast<PetscObject>(snes);
56
57 PetscReal residual_min, residual_max;
58
1/2
✓ Branch 2 taken 20911 times.
✗ Branch 3 not taken.
 20911 p_linear_problem->vector().min(&residual_min);
59
1/2
✓ Branch 2 taken 20911 times.
✗ Branch 3 not taken.
 20911 p_linear_problem->vector().max(&residual_max);
60
61 int tablevel;
62
1/2
✓ Branch 1 taken 20911 times.
✗ Branch 2 not taken.
 20911 PetscObjectGetTabLevel(snes_object,&tablevel);
63
2/4
✓ Branch 1 taken 20911 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20911 times.
✗ Branch 5 not taken.
 20911 PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_WORLD, tablevel);
64
2/4
✓ Branch 1 taken 20911 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20911 times.
✗ Branch 5 not taken.
 20911 PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,
65 "%d SNES Function Norm is %14.12e and displacement extrema are %8.6e and %8.6e. Residual between %8.6e and %8.6e\n",
66 its,
67 norm,
68 evaluation_state_Min, evaluation_state_Max,
69 residual_min, residual_max
70 );
71
72
2/4
✓ Branch 1 taken 20911 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20911 times.
✗ Branch 5 not taken.
 20911 PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_WORLD, tablevel);
73
74 20911 PetscFunctionReturn(0);
75 }
76 } // namespace anonymous
77
78 /***********************************************************************************/
79 /***********************************************************************************/
80
81 518 SNESInterface::SNESInterface()
82 {
83 518 }
84
85 /***********************************************************************************/
86 /***********************************************************************************/
87
88 518 SNESInterface::~SNESInterface()
89 {
90
2/2
✓ Branch 0 taken 283 times.
✓ Branch 1 taken 235 times.
 518 if (m_doUseSNES) {
91 283 SNESDestroy(&m_SNES);
92 }
93 518 }
94
95 /***********************************************************************************/
96 /***********************************************************************************/
97
98 283 void SNESInterface::init()
99 {
100 283 init( MpiInfo::petscComm() );
101 283 }
102
103 /***********************************************************************************/
104 /***********************************************************************************/
105
106 283 void SNESInterface::init(MPI_Comm comm)
107 {
108 283 m_comm = comm;
109
110 // First case of Newton-Raphson method such as implemented by Petsc.
111 // By extracting the KSP and PC contexts from the SNES context, we can then
112 // directly call any KSP, and PC routines to set various options.
113
1/2
✓ Branch 1 taken 283 times.
✗ Branch 2 not taken.
 283 SNESCreate(m_comm, &m_SNES);
114
115 // We deactivate line-search in order to use standard Newton-Raphson
116
1/2
✓ Branch 1 taken 283 times.
✗ Branch 2 not taken.
 283 SNESSetType(m_SNES,SNESNEWTONLS);
117 SNESLineSearch ls;
118
1/2
✓ Branch 1 taken 283 times.
✗ Branch 2 not taken.
 283 SNESGetLineSearch(m_SNES, &ls);
119
1/2
✓ Branch 1 taken 283 times.
✗ Branch 2 not taken.
 283 SNESLineSearchSetType(ls,SNESLINESEARCHBASIC);
120 283 }
121
122 /***********************************************************************************/
123 /***********************************************************************************/
124
125 7242 void SNESInterface::getSolutionUpdate(PetscVector& rDeltaX)
126 {
127 Vec increment;
128 // Warning: SNESGetSolutionUpdate returns -delta U
129
1/2
✓ Branch 1 taken 7242 times.
✗ Branch 2 not taken.
 7242 SNESGetSolutionUpdate(m_SNES, &increment);
130
131 // Copy to our vector
132
2/4
✓ Branch 1 taken 7242 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7242 times.
✗ Branch 5 not taken.
 7242 VecCopy(increment,rDeltaX.toPetsc());
133 7242 }
134
135 /***********************************************************************************/
136 /***********************************************************************************/
137
138 247 void SNESInterface::setFunction(PetscVector& r, PetscErrorCode (*f)(SNES,Vec,Vec,void*), void* ctx)
139 {
140 247 SNESSetFunction(m_SNES, r.toPetsc(), f, ctx);
141 247 }
142
143 /***********************************************************************************/
144 /***********************************************************************************/
145
146 247 void SNESInterface::setJacobian(
147 PetscMatrix& Amat,
148 PetscMatrix& Pmat,
149 PetscErrorCode (*J)(SNES,Vec,Mat,Mat,void*),
150 void* ctx
151 )
152 {
153 247 SNESSetJacobian(m_SNES, Amat.toPetsc(), Pmat.toPetsc(), J, ctx);
154 247 }
155
156 /***********************************************************************************/
157 /***********************************************************************************/
158
159 283 void SNESInterface::getKSP(KSP* pKSP)
160 {
161 283 SNESGetKSP(m_SNES, pKSP);
162 283 }
163
164 /***********************************************************************************/
165 /***********************************************************************************/
166
167 283 void SNESInterface::setMonitorSetLinearProblem(LinearProblem* pLinearProblem)
168 {
169 283 SNESMonitorSet(m_SNES, snesViewerLinearProblem, pLinearProblem, FELISCE_PETSC_NULLPTR);
170 283 }
171
172 /***********************************************************************************/
173 /***********************************************************************************/
174
175 283 void SNESInterface::setTolerances(
176 const double absoluteTolerance,
177 const double relativeTolerance,
178 const double solutionTolerance,
179 const int maxIterationSNES,
180 const int maxFunctionEvaluatedSNES,
181 const double divergenceTolerance
182 )
183 {
184 283 SNESSetTolerances(m_SNES,
185 absoluteTolerance,
186 relativeTolerance,
187 solutionTolerance,
188 maxIterationSNES,
189 maxFunctionEvaluatedSNES
190 );
191
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 279 times.
 283 if (divergenceTolerance > 0.0) {
192 4 SNESSetDivergenceTolerance(m_SNES, divergenceTolerance);
193 }
194 283 }
195
196 /***********************************************************************************/
197 /***********************************************************************************/
198
199 7095 PetscInt SNESInterface::solve(PetscVector& rX)
200 {
201
2/4
✓ Branch 1 taken 7095 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7095 times.
✗ Branch 5 not taken.
 7095 SNESSolve(m_SNES, FELISCE_PETSC_NULLPTR, rX.toPetsc());
202
203 SNESConvergedReason cvg_reason;
204
1/2
✓ Branch 1 taken 7095 times.
✗ Branch 2 not taken.
 7095 SNESGetConvergedReason(m_SNES,&cvg_reason);
205 PetscInt its;
206
1/2
✓ Branch 1 taken 7095 times.
✗ Branch 2 not taken.
 7095 SNESGetIterationNumber(m_SNES,&its);
207
3/4
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 7063 times.
✓ Branch 2 taken 32 times.
✗ Branch 3 not taken.
 7095 if (cvg_reason > 0 || cvg_reason == SNES_DIVERGED_MAX_IT) { // NOTE: See https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/SNES/SNESConvergedReason.html
208
4/6
✓ Branch 0 taken 2226 times.
✓ Branch 1 taken 4293 times.
✓ Branch 2 taken 544 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 32 times.
✗ Branch 5 not taken.
 7095 switch (cvg_reason) {
209
1/2
✓ Branch 1 taken 2226 times.
✗ Branch 2 not taken.
 2226 case SNES_CONVERGED_FNORM_ABS: PetscPrintf(m_comm, "SNES converged in %d iterations due to ||F|| < atol \n", its); break;
210
1/2
✓ Branch 1 taken 4293 times.
✗ Branch 2 not taken.
 4293 case SNES_CONVERGED_FNORM_RELATIVE: PetscPrintf(m_comm, "SNES converged in %d iterations due to ||F|| < rtol*||F_initial|| \n", its); break;
211
1/2
✓ Branch 1 taken 544 times.
✗ Branch 2 not taken.
 544 case SNES_CONVERGED_SNORM_RELATIVE: PetscPrintf(m_comm, "SNES converged in %d iterations due to Newton computed step size small; || delta x || < stol || x || \n", its); break;
212 ✗ case SNES_CONVERGED_ITS: PetscPrintf(m_comm, "SNES converged in %d iterations due to maximum iterations reached \n", its); break;
213
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 32 case SNES_DIVERGED_MAX_IT: PetscPrintf(m_comm, "SNES diverged due to SNESSolve() has reached the maximum number of iterations requested: %d. Calculations continues as normal \n", its); break;
214 ✗ default: PetscPrintf(m_comm, "SNES converged in %d iterations \n", its);
215 }
216 } else {
217 ✗ std::string error_message = "Convergence failed in SNES solver. ";
218 ✗ switch (cvg_reason) {
219 // case SNES_CONVERGED_ITERATING: error_message += "SNES diverged due to maximum iterations reached"; break;
220 ✗ case SNES_DIVERGED_FUNCTION_DOMAIN: error_message += "SNES diverged due to the new x location passed the function is not in the domain of F"; break;
221 ✗ case SNES_DIVERGED_FUNCTION_COUNT: error_message += "SNES diverged due to the user provided function has been called more times then the final argument to SNESSetTolerances() "; break;
222 ✗ case SNES_DIVERGED_LINEAR_SOLVE: error_message += "SNES diverged due to the linear solve failed"; break;
223 ✗ case SNES_DIVERGED_FNORM_NAN: error_message += "SNES diverged due to the 2-norm of the current function evaluation is not-a-number (NaN), this is usually caused by a division of 0 by 0"; break;
224 // case SNES_DIVERGED_MAX_IT: error_message += "SNES diverged due to SNESSolve() has reached the maximum number of iterations requested"; break;
225 ✗ case SNES_DIVERGED_LINE_SEARCH: error_message += "SNES diverged due to the line search failed "; break;
226 ✗ case SNES_DIVERGED_INNER: error_message += "SNES diverged due to inner solve failed"; break;
227 ✗ case SNES_DIVERGED_LOCAL_MIN: error_message += "SNES diverged due to || J^T b || is small, implies converged to local minimum of F()"; break;
228 ✗ case SNES_DIVERGED_DTOL: error_message += "SNES diverged due to || F || > divtol*||F_initial||"; break;
229 ✗ case SNES_DIVERGED_JACOBIAN_DOMAIN: error_message += "SNES diverged due to Jacobian calculation does not make sense"; break;
230 ✗ case SNES_DIVERGED_TR_DELTA: error_message += "SNES diverged due to SNES_DIVERGED_TR_DELTA"; break;
231 ✗ default: error_message += "";
232 }
233 ✗ FEL_ERROR(error_message);
234 }
235
236 7095 return its;
237 }
238
239 }
240