Skip to content

Commit 31539bb

Browse files
yib0liuyib0liu
committed
Add Trilinos spec
1 parent 8f0210a commit 31539bb

File tree

5 files changed

+67
-165
lines changed

5 files changed

+67
-165
lines changed

‎linear-solver-spec.json‎

Lines changed: 27 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -160,8 +160,9 @@
160160
"default": null,
161161
"type": "object",
162162
"optional": [
163-
"solver",
164-
"precond"
163+
"max_iter",
164+
"tolerance",
165+
"block_size"
165166
],
166167
"doc": "Settings for the Trilinos solver."
167168
},
@@ -433,5 +434,29 @@
433434
"default": 0,
434435
"type": "float",
435436
"doc": "Aggregation epsilon strong."
437+
},
438+
{
439+
"pointer": "/Trilinos/max_iter",
440+
"default": 1000,
441+
"type": "int",
442+
"doc": "Maximum number of iterations."
443+
},
444+
{
445+
"pointer": "/Trilinos/block_size",
446+
"default": 3,
447+
"type": "int",
448+
"doc": "Aggregation epsilon strong."
449+
},
450+
{
451+
"pointer": "/Trilinos/tolerance",
452+
"default": 1e-8,
453+
"type": "float",
454+
"doc": "Convergence tolerance."
455+
},
456+
{
457+
"pointer": "/Trilinos/is_nullspace",
458+
"default": false,
459+
"type": "bool",
460+
"doc": "Is nullspace or not."
436461
}
437462
]

‎src/polysolve/linear/Solver.cpp‎

Lines changed: 0 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -45,33 +45,6 @@
4545
#include <unsupported/Eigen/IterativeSolvers>
4646

4747
////////////////////////////////////////////////////////////////////////////////
48-
Eigen::MatrixXd test_vertices;
49-
Eigen::MatrixXd init_vertices;
50-
std::vector<int> test_boundary_nodes;
51-
Eigen::MatrixXd remove_boundary_vertices(const Eigen::MatrixXd &vertices, const std::vector<int> &boundary_nodes)
52-
{
53-
// Remove boundary vertices
54-
if (boundary_nodes.empty())
55-
{
56-
return vertices;
57-
}
58-
else
59-
{
60-
std::vector<int> order_nodes = boundary_nodes;
61-
std::sort(order_nodes.begin(), order_nodes.end());
62-
Eigen::MatrixXd out_vertices;
63-
std::vector<int> keep;
64-
for (int i = 0; i < vertices.rows(); i++)
65-
{
66-
if (!std::binary_search(order_nodes.begin(), order_nodes.end(),i))
67-
{
68-
keep.push_back(i);
69-
}
70-
}
71-
out_vertices = vertices(keep, Eigen::all);
72-
return out_vertices;
73-
}
74-
}
7548

7649
namespace polysolve::linear
7750
{

‎src/polysolve/linear/Solver.hpp‎

Lines changed: 0 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -14,11 +14,6 @@
1414
Base(const Base &) = delete; \
1515
Base &operator=(const Base &) = delete;
1616

17-
extern Eigen::MatrixXd init_vertices;
18-
extern Eigen::MatrixXd test_vertices;
19-
extern std::vector<int> test_boundary_nodes;
20-
Eigen::MatrixXd remove_boundary_vertices(const Eigen::MatrixXd &vertices, const std::vector<int> &boundary_nodes);
21-
2217
////////////////////////////////////////////////////////////////////////////////
2318
// TODO:
2419
// - [ ] Support both RowMajor + ColumnMajor sparse matrices
@@ -34,11 +29,6 @@ namespace spdlog
3429

3530
namespace polysolve::linear
3631
{
37-
#ifdef POLYSOLVE_LARGE_INDEX
38-
typedef Eigen::SparseMatrix<double, Eigen::ColMajor, std::ptrdiff_t> StiffnessMatrix;
39-
#else
40-
typedef Eigen::SparseMatrix<double, Eigen::ColMajor> StiffnessMatrix;
41-
#endif
4232
/**
4333
* @brief Base class for linear solver.
4434
*/

‎src/polysolve/linear/TrilinosSolver.cpp‎

Lines changed: 35 additions & 126 deletions
Original file line numberDiff line numberDiff line change
@@ -5,103 +5,42 @@
55
#include <string>
66
#include <vector>
77
#include <unsupported/Eigen/SparseExtra>
8-
/////////////////////////////////s///////////////////////////////////////////////
8+
////////////////////////////////////////////////////////////////////////////////
99

10-
namespace polysolve::linear
10+
Eigen::MatrixXd test_vertices;
11+
Eigen::MatrixXd init_vertices;
12+
std::vector<int> test_boundary_nodes;
13+
Eigen::MatrixXd remove_boundary_vertices(const Eigen::MatrixXd &vertices, const std::vector<int> &boundary_nodes)
1114
{
12-
namespace{
13-
////////////////////////////////////////////////////////////////
14-
// int rigid_body_mode(int ndim, const std::vector<double> &coo, std::vector<double> &B, bool transpose = true) {
15-
16-
// size_t n = coo.size();
17-
// int nmodes = (ndim == 2 ? 3 : 6);
18-
// B.resize(n * nmodes, 0.0);
19-
20-
// const int stride1 = transpose ? 1 : nmodes;
21-
// const int stride2 = transpose ? n : 1;
22-
// // int stride1 = nmodes;
23-
// // int stride2 = 1;
24-
25-
// double sn = 1 / sqrt(n);
26-
27-
// if (ndim == 2) {
28-
// for(size_t i = 0; i < n; ++i) {
29-
// size_t nod = i / ndim;
30-
// size_t dim = i % ndim;
31-
32-
// double x = coo[nod * 2 + 0];
33-
// double y = coo[nod * 2 + 1];
34-
35-
// // Translation
36-
// B[i * stride1 + dim * stride2] = sn;
37-
38-
// // Rotation
39-
// switch(dim) {
40-
// case 0:
41-
// B[i * stride1 + 2 * stride2] = -y;
42-
// break;
43-
// case 1:
44-
// B[i * stride1 + 2 * stride2] = x;
45-
// break;
46-
// }
47-
// }
48-
// } else if (ndim == 3) {
49-
// for(size_t i = 0; i < n; ++i) {
50-
// size_t nod = i / ndim;
51-
// size_t dim = i % ndim;
52-
53-
// double x = coo[nod * 3 + 0];
54-
// double y = coo[nod * 3 + 1];
55-
// double z = coo[nod * 3 + 2];
56-
57-
// // Translation
58-
// B[i * stride1 + dim * stride2] = sn;
59-
60-
// // Rotation
61-
// switch(dim) {
62-
// case 0:
63-
// B[i * stride1 + 5 * stride2] = -y;
64-
// B[i * stride1 + 4 * stride2] = z;
65-
// break;
66-
// case 1:
67-
// B[i * stride1 + 5 * stride2] = x;
68-
// B[i * stride1 + 3 * stride2] = -z;
69-
// break;
70-
// case 2:
71-
// B[i * stride1 + 3 * stride2] = y;
72-
// B[i * stride1 + 4 * stride2] = -x;
73-
// break;
74-
// }
75-
// }
76-
// }
77-
78-
// // Orthonormalization
79-
// std::array<double, 6> dot;
80-
// for(int i = ndim; i < nmodes; ++i) {
81-
// std::fill(dot.begin(), dot.end(), 0.0);
82-
// for(size_t j = 0; j < n; ++j) {
83-
// for(int k = 0; k < i; ++k)
84-
// dot[k] += B[j * stride1 + k * stride2] * B[j * stride1 + i * stride2];
85-
// }
86-
// double s = 0.0;
87-
// for(size_t j = 0; j < n; ++j) {
88-
// for(int k = 0; k < i; ++k)
89-
// B[j * stride1 + i * stride2] -= dot[k] * B[j * stride1 + k * stride2];
90-
// s += B[j * stride1 + i * stride2] * B[j * stride1 + i * stride2];
91-
// }
92-
// s = sqrt(s);
93-
// for(size_t j = 0; j < n; ++j)
94-
// B[j * stride1 + i * stride2] /= s;
95-
// }
96-
// return nmodes;
97-
// }
98-
// Produce a vector of rigid body modes
99-
15+
// Remove boundary vertices
16+
if (boundary_nodes.empty())
17+
{
18+
return vertices;
10019
}
20+
else
21+
{
22+
std::vector<int> order_nodes = boundary_nodes;
23+
std::sort(order_nodes.begin(), order_nodes.end());
24+
Eigen::MatrixXd out_vertices;
25+
std::vector<int> keep;
26+
for (int i = 0; i < vertices.rows(); i++)
27+
{
28+
if (!std::binary_search(order_nodes.begin(), order_nodes.end(),i))
29+
{
30+
keep.push_back(i);
31+
}
32+
}
33+
out_vertices = vertices(keep, Eigen::all);
34+
return out_vertices;
35+
}
36+
}
37+
38+
namespace polysolve::linear
39+
{
10140
TrilinosSolver::TrilinosSolver()
10241
{
10342
precond_num_ = 0;
104-
#ifdef HAVE_MPI
43+
// #ifdef HAVE_MPI
10544
int done_already;
10645
MPI_Initialized(&done_already);
10746
if (!done_already)
@@ -114,9 +53,9 @@ namespace polysolve::linear
11453
MPI_Init(&argc, &argvv);
11554
CommPtr = new Epetra_MpiComm(MPI_COMM_WORLD);
11655
}
117-
#else
118-
CommPtr=new Epetra_SerialComm;
119-
#endif
56+
// #else
57+
// CommPtr=new Epetra_SerialComm;
58+
// #endif
12059
}
12160

12261
////////////////////////////////////////////////////////////////
@@ -226,36 +165,6 @@ namespace polysolve::linear
226165
Teuchos::ParameterList MLList;
227166
TrilinosML_SetDefaultOptions(MLList);
228167
MLList.set("PDE equations",numPDEs);
229-
230-
//Set null space
231-
232-
// if (true)
233-
// {
234-
// int n=test_vertices.rows();
235-
// int NRbm=0;
236-
// int NscalarDof=0;
237-
238-
// if (numPDEs==2)
239-
// {
240-
// NRbm=3;
241-
// rbm=new double[n*(NRbm+NscalarDof)*(numPDEs+NscalarDof)];
242-
// std::vector<double> z_coord(n,0);
243-
// ML_Coord2RBM(n,test_vertices.col(0).data(),test_vertices.col(1).data(),z_coord.data(),rbm,numPDEs,NscalarDof);
244-
// }
245-
// else
246-
// {
247-
// NRbm=6;
248-
// rbm=new double[n*(NRbm+NscalarDof)*(numPDEs+NscalarDof)];
249-
// ML_Coord2RBM(n,test_vertices.col(0).data(),test_vertices.col(1).data(),test_vertices.col(2).data(),rbm,numPDEs,NscalarDof);
250-
// }
251-
252-
// MLList.set("null space: vectors",rbm);
253-
// MLList.set("null space: dimension", NRbm);
254-
// MLList.set("null space: type", "pre-computed");
255-
// MLList.set("aggregation: threshold",0.00);
256-
// }
257-
258-
///////////////////////////////////////////////////////////////////////
259168

260169
if (is_nullspace_)
261170
{
@@ -342,9 +251,9 @@ namespace polysolve::linear
342251
delete A;
343252
delete rowMap;
344253
delete MLPrec;
345-
#ifdef HAVE_MPI
254+
// #ifdef HAVE_MPI
346255
MPI_Finalize() ;
347-
#endif
256+
// #endif
348257
}
349258
}
350259

‎src/polysolve/linear/TrilinosSolver.hpp‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,11 @@
4141
// column-major matrix, the solver will actually solve A^T x = b.
4242
//
4343

44+
extern Eigen::MatrixXd init_vertices;
45+
extern Eigen::MatrixXd test_vertices;
46+
extern std::vector<int> test_boundary_nodes;
47+
Eigen::MatrixXd remove_boundary_vertices(const Eigen::MatrixXd &vertices, const std::vector<int> &boundary_nodes);
48+
4449
namespace polysolve::linear
4550
{
4651

0 commit comments

Comments
 (0)