GCC Code Coverage Report

Directory: ./
File: Tools/ksp_solver.hpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 103 118 87.3%
Branches: 101 314 32.2%
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 #include <petscmat.h>
19 #include <petscksp.h>
20 #include "Core/NoThirdPartyWarning/Petsc/sys.hpp" // Jacques 13/07/22: added due to definition of FELISCE_PETSC_NULLPTR
21
22 // Project includes
23 #include "Core/filesystemUtil.hpp"
24
25
26 #if !defined(FELISCE_KSP_SOLVER_H_INCLUDED )
27 #define FELISCE_KSP_SOLVER_H_INCLUDED
28
29 72 int SolveMatrix(
30 const std::string FilenameMatrix,
31 const std::string FilenameVector,
32 const std::string SolverType,
33 const std::string PreconditionerType = "none"
34 )
35 {
36 PetscErrorCode ierr;
37 PetscInt its,m,n,mvec;
38 PetscReal norm;
39 Vec x,b,u;
40 Mat A;
41 KSP ksp;
42 PetscViewer fd;
43 PetscLogDouble t1, t2;
44 72 const PetscReal rtol = 1.e-8;
45
46
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = PetscInitialize(FELISCE_PETSC_NULLPTR, FELISCE_PETSC_NULLPTR, FELISCE_PETSC_NULLPTR, FELISCE_PETSC_NULLPTR );
47
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
48
49 /* Read matrix and RHS */
50 72 const char* file_matrix = FilenameMatrix.c_str();
51
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file_matrix,FILE_MODE_READ,&fd);
52
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
53
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = MatCreate(PETSC_COMM_WORLD,&A);
54
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
55
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = MatSetType(A,MATMPIAIJ);
56
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
57
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = MatLoad(A,fd);
58
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
59
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = PetscViewerDestroy(&fd);
60
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
61
62 72 const char* file_vector = FilenameVector.c_str();
63
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file_vector,FILE_MODE_READ,&fd);
64
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
65
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = VecCreate(PETSC_COMM_WORLD,&b);
66
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
67
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = VecSetType(b,VECMPI);
68
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
69
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = VecLoad(b,fd);
70
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
71
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = PetscViewerDestroy(&fd);
72
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
73
74 /*
75 If the load matrix is larger then the vector, due to being padded
76 to match the blocksize then create a new padded vector
77 */
78
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = MatGetSize(A,&m,&n);
79
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
80
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = VecGetSize(b,&mvec);
81
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
82
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
 72 if (m > mvec) {
83 Vec tmp;
84 PetscScalar *bold,*bnew;
85 /* create a new vector b by padding the old one */
86 ✗ ierr = VecCreate(PETSC_COMM_WORLD,&tmp);
87 ✗ CHKERRQ(ierr);
88 ✗ ierr = VecSetSizes(tmp,PETSC_DECIDE,m);
89 ✗ CHKERRQ(ierr);
90 ✗ ierr = VecSetFromOptions(tmp);
91 ✗ CHKERRQ(ierr);
92 ✗ ierr = VecGetArray(tmp,&bnew);
93 ✗ CHKERRQ(ierr);
94 ✗ ierr = VecGetArray(b,&bold);
95 ✗ CHKERRQ(ierr);
96 ✗ ierr = PetscArraycpy(bnew,bold,mvec);
97 ✗ CHKERRQ(ierr);
98 ✗ ierr = VecDestroy(&b);
99 ✗ CHKERRQ(ierr);
100 ✗ b = tmp;
101 }
102
103 /* Set up solution */
104
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = VecDuplicate(b,&x);
105
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
106
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = VecDuplicate(b,&u);
107
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
108
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = VecSet(x,0.0);
109
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
110
111 /* Solve system */
112
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);
113
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
114
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = KSPSetOperators(ksp,A,A);
115
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
116
117 // Retrieve the preconditioner
118 PC pc;
119
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 KSPGetPC(ksp, &pc);
120
121 // Assigning the solver
122
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 68 times.
 72 if (SolverType == "MUMPS") {
123
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 KSPSetType(ksp, KSPPREONLY);
124
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 PCSetType(pc, PCLU);
125
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 PCFactorSetMatSolverType(pc, MATSOLVERMUMPS);
126 // } else if (SolverType == "PASTIX") { // NOTE: Does not work properly yet with PETSc integration
127 // KSPSetType(ksp, KSPPREONLY);
128 // PC pc;
129 // KSPGetPC(ksp, &pc);
130 // PCSetType(pc, PCLU);
131 // PCFactorSetMatSolverType(pc, MATSOLVERPASTIX);
132 #if SUPERLU_DIST_COMPILED
133
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 64 times.
 68 } else if (SolverType == "SUPERLU_DIST") {
134
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 KSPSetType(ksp, KSPPREONLY);
135 PC pc;
136
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 KSPGetPC(ksp, &pc);
137
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 PCSetType(pc, PCLU);
138
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 PCFactorSetMatSolverType(pc, MATSOLVERSUPERLU_DIST);
139 #endif
140 #if HYPRE_COMPILED
141
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 60 times.
 64 } else if (SolverType == "hypre") {
142
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 KSPSetType(ksp, "gmres");
143
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 PCSetType(pc, "hypre");
144
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 PCHYPRESetType(pc, PreconditionerType.c_str());
145 #endif
146 } else {
147
1/2
✓ Branch 2 taken 60 times.
✗ Branch 3 not taken.
 60 KSPSetType(ksp, SolverType.c_str());
148
1/2
✓ Branch 2 taken 60 times.
✗ Branch 3 not taken.
 60 PCSetType(pc, PreconditionerType.c_str());
149 }
150
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 PetscOptionsSetValue(NULL, "-mat_mumps_icntl_20", "0");
151
152
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
153
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = KSPSetTolerances(ksp, rtol, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
154
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
155
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = KSPSetFromOptions(ksp);
156
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
157
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = PetscTime(&t1);
158
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
159
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = KSPSolve(ksp,b,x);
160
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
161
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = PetscTime(&t2);
162
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
163
164 /* Show result */
165
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = MatMult(A,x,u);
166
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
167
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = VecAXPY(u,-1.0,b);
168
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
169
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = VecNorm(u,NORM_2,&norm);
170
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
171
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = KSPGetIterationNumber(ksp,&its);
172
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
173 72 const char* solver_type = SolverType.c_str();
174
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = PetscPrintf(PETSC_COMM_WORLD, "Solver = %s \n", solver_type);
175
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = PetscPrintf(PETSC_COMM_WORLD, "Elapse time: %lf s \n", t2 - t1);
176
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
177 KSPConvergedReason reason;
178
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = KSPGetConvergedReason(ksp, &reason);
179
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
180
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = PetscPrintf(PETSC_COMM_WORLD, "Converged reason: %s \n", KSPConvergedReasons[reason]);
181
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
182
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = PetscPrintf(PETSC_COMM_WORLD, "Number of iterations = %3d \n",its);
183
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
184
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = PetscPrintf(PETSC_COMM_WORLD, "Residual norm %g \n",(double)norm);
185
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
186
2/4
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 72 times.
✗ Branch 5 not taken.
 72 ierr = KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD);
187
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
188
189 /* Cleanup */
190
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = KSPDestroy(&ksp);
191
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
192
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = VecDestroy(&x);
193
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
194
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = VecDestroy(&b);
195
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
196
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = VecDestroy(&u);
197
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
198
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = MatDestroy(&A);
199
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 72 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 72 CHKERRQ(ierr);
200
201
1/2
✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
 72 ierr = PetscFinalize();
202 72 return ierr;
203 }
204
205 #endif
206