Skip to content

Commit 5650bc9

Browse files
akva2akva2
committed
added: support for K-refinement cycles
1 parent f815348 commit 5650bc9

File tree

14 files changed

+648
-34
lines changed

14 files changed

+648
-34
lines changed

‎Apps/Common/SIMKCyclePC.h‎

Lines changed: 228 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,228 @@
1+
//==============================================================================
2+
//!
3+
//! \file SIMKCyclePC.h
4+
//!
5+
//! \date Nov 28 2012
6+
//!
7+
//! \author Arne Morten Kvarving / SINTEF
8+
//!
9+
//! \brief Use a simulator as a preconditioner for K-refinement cycles.
10+
//!
11+
//==============================================================================
12+
#ifndef SIM_K_CYCLE_PC_H_
13+
#define SIM_K_CYCLE_PC_H_
14+
15+
#ifdef HAS_PETSC
16+
17+
#include "ASMs2D.h"
18+
#include "ASMs3D.h"
19+
#include "SAM.h"
20+
#include "PETScSolParams.h"
21+
#include <GoTools/geometry/SplineSurface.h>
22+
#include <GoTools/trivariate/SplineVolume.h>
23+
24+
25+
/*!
26+
\brief Template for K-cycle preconditioning.
27+
\details Currently only single patch models.
28+
*/
29+
30+
template<class Sim>
31+
class SIMKCyclePC : public PETScPC
32+
{
33+
public:
34+
//! \brief Default constructor.
35+
//! \param fromSim The model for the solution
36+
//! \param pcSim The model for the preconditioner
37+
SIMKCyclePC(Sim& fromSim, Sim& pcSim) :
38+
S1(pcSim), fSim(fromSim),
39+
B(SparseMatrix::SUPERLU), BT(SparseMatrix::SUPERLU)
40+
{
41+
setup();
42+
}
43+
44+
//! \brief Setup the restriction/prolongation operators.
45+
//! \details Operators are constructed through Greville point projection.
46+
void setup()
47+
{
48+
if (fSim.getPatch(1)->getNoSpaceDim() == 2) {
49+
const ASMs2D* fpch = dynamic_cast<const ASMs2D*>(fSim.getPatch(1));
50+
const ASMs2D* pch = dynamic_cast<const ASMs2D*>(S1.getPatch(1));
51+
std::array<RealArray,2> gPrm;
52+
fpch->getGrevilleParameters(gPrm[0], 0);
53+
fpch->getGrevilleParameters(gPrm[1], 1);
54+
N.resize(gPrm[0].size()*gPrm[1].size(), S1.getNoDOFs());
55+
NT.resize(S1.getNoDOFs(), gPrm[0].size()*gPrm[1].size());
56+
std::vector<Go::BasisDerivsSf> spline;
57+
pch->getBasis(1)->computeBasisGrid(gPrm[0], gPrm[1], spline);
58+
int p1 = pch->getBasis(1)->order_u();
59+
int p2 = pch->getBasis(1)->order_v();
60+
int n1 = pch->getBasis(1)->numCoefs_u();
61+
int n2 = pch->getBasis(1)->numCoefs_v();
62+
for (size_t ip = 0; ip < gPrm[0].size()*gPrm[1].size(); ++ip) {
63+
IntVec idx;
64+
ASMs2D::scatterInd(n1,n2,p1,p2,spline[ip].left_idx,idx);
65+
for (size_t k = 0; k < spline[ip].basisValues.size(); ++k) {
66+
N(ip+1, idx[k]+1) += spline[ip].basisValues[k];
67+
NT(idx[k]+1, ip+1) += spline[ip].basisValues[k];
68+
}
69+
}
70+
fpch->getBasis(1)->computeBasisGrid(gPrm[0], gPrm[1], spline);
71+
p1 = fpch->getBasis(1)->order_u();
72+
p2 = fpch->getBasis(1)->order_v();
73+
n1 = fpch->getBasis(1)->numCoefs_u();
74+
n2 = fpch->getBasis(1)->numCoefs_v();
75+
B.resize(fSim.getNoDOFs(), fSim.getNoDOFs());
76+
BT.resize(fSim.getNoDOFs(), fSim.getNoDOFs());
77+
78+
for (size_t ip = 0; ip < gPrm[0].size()*gPrm[1].size(); ++ip) {
79+
IntVec idx;
80+
ASMs2D::scatterInd(n1,n2,p1,p2,spline[ip].left_idx,idx);
81+
for (size_t k = 0; k < spline[ip].basisValues.size(); ++k) {
82+
B(ip+1, idx[k]+1) += spline[ip].basisValues[k];
83+
BT(idx[k]+1, ip+1) += spline[ip].basisValues[k];
84+
}
85+
}
86+
} else {
87+
const ASMs3D* fpch = dynamic_cast<const ASMs3D*>(fSim.getPatch(1));
88+
const ASMs3D* pch = dynamic_cast<const ASMs3D*>(S1.getPatch(1));
89+
std::array<RealArray,3> gPrm;
90+
fpch->getGrevilleParameters(gPrm[0], 0);
91+
fpch->getGrevilleParameters(gPrm[1], 1);
92+
fpch->getGrevilleParameters(gPrm[2], 2);
93+
N.resize(gPrm[0].size()*gPrm[1].size()*gPrm[2].size(), S1.getNoDOFs());
94+
NT.resize(S1.getNoDOFs(), gPrm[0].size()*gPrm[1].size()*gPrm[2].size());
95+
std::vector<Go::BasisDerivs> spline;
96+
pch->getBasis(1)->computeBasisGrid(gPrm[0], gPrm[1], gPrm[2], spline);
97+
int p1 = pch->getBasis(1)->order(0);
98+
int p2 = pch->getBasis(1)->order(1);
99+
int p3 = pch->getBasis(1)->order(2);
100+
int n1 = pch->getBasis(1)->numCoefs(0);
101+
int n2 = pch->getBasis(1)->numCoefs(1);
102+
int n3 = pch->getBasis(1)->numCoefs(2);
103+
for (size_t ip = 0; ip < gPrm[0].size()*gPrm[1].size()*gPrm[2].size(); ++ip) {
104+
IntVec idx;
105+
ASMs3D::scatterInd(n1,n2,n3,p1,p2,p3,spline[ip].left_idx,idx);
106+
for (size_t k = 0; k < spline[ip].basisValues.size(); ++k) {
107+
N(ip+1, idx[k]+1) += spline[ip].basisValues[k];
108+
NT(idx[k]+1, ip+1) += spline[ip].basisValues[k];
109+
}
110+
}
111+
fpch->getBasis(1)->computeBasisGrid(gPrm[0], gPrm[1], gPrm[2], spline);
112+
p1 = fpch->getBasis(1)->order(0);
113+
p2 = fpch->getBasis(1)->order(1);
114+
p3 = fpch->getBasis(1)->order(2);
115+
n1 = fpch->getBasis(1)->numCoefs(0);
116+
n2 = fpch->getBasis(1)->numCoefs(1);
117+
n3 = fpch->getBasis(1)->numCoefs(2);
118+
B.resize(fSim.getNoDOFs(), fSim.getNoDOFs());
119+
BT.resize(fSim.getNoDOFs(), fSim.getNoDOFs());
120+
121+
for (size_t ip = 0; ip < gPrm[0].size()*gPrm[1].size()*gPrm[2].size(); ++ip) {
122+
IntVec idx;
123+
ASMs3D::scatterInd(n1,n2,n3,p1,p2,p3,spline[ip].left_idx,idx);
124+
for (size_t k = 0; k < spline[ip].basisValues.size(); ++k) {
125+
B(ip+1, idx[k]+1) += spline[ip].basisValues[k];
126+
BT(idx[k]+1, ip+1) += spline[ip].basisValues[k];
127+
}
128+
}
129+
}
130+
}
131+
132+
//! \brief Evaluate preconditioner.
133+
bool eval(Vec& x, Vec& y) override
134+
{
135+
// copy vector into a standard vector for expansion
136+
int len;
137+
VecGetSize(x, &len);
138+
StdVector X(len);
139+
std::vector<PetscInt> idx(len);
140+
std::iota(idx.begin(), idx.end(), 0);
141+
VecGetValues(x, len, idx.data(), X.getPtr());
142+
143+
// Expand solution vector x to DOF vector
144+
StdVector dofSol(fSim.getNoDOFs());
145+
fSim.getSAM()->expandSolution(X, dofSol);
146+
147+
// Evaluate restriction - N^T * (B^T)^-1 * dofSol
148+
BT.solve(dofSol);
149+
StdVector restrict(NT.rows());
150+
NT.multiply(dofSol, restrict);
151+
152+
// copy to solution vector
153+
StdVector* eqsVec = S1.getVector();
154+
eqsVec->init();
155+
vectorToSolVec(S1, *eqsVec, restrict);
156+
eqsVec->beginAssembly();
157+
eqsVec->endAssembly();
158+
159+
// Perform inner solve, store in restrict
160+
double cond;
161+
int oldLev = this->S1.msgLevel;
162+
this->S1.msgLevel = 0;
163+
if (!this->S1.solveSystem(restrict, 0, &cond, "displacement", false))
164+
return false;
165+
this->S1.msgLevel = oldLev;
166+
167+
// Evaluate prolonged solution - dofSol = B^-1 * N * restrict
168+
N.multiply(restrict, dofSol);
169+
B.solve(dofSol);
170+
171+
// Inject into solution vector y
172+
const int* meqn = fSim.getSAM()->getMEQN();
173+
for (size_t i = 1; i <= fSim.getPatch(1)->getNoNodes(); ++i) {
174+
int eq = meqn[i-1];
175+
if (eq > 0)
176+
VecSetValue(y, eq-1, dofSol(i), INSERT_VALUES);
177+
}
178+
179+
return true;
180+
}
181+
182+
//! \brief Evaluate the right-hand-side system - the prolongiation of the coarse system solution.
183+
//! \param rhs The resulting vector
184+
void evalRHS(PETScVector& rhs)
185+
{
186+
// Solve coarse system, store in solCoarse
187+
StdVector solCoarse(S1.getNoDOFs());
188+
S1.solveSystem(solCoarse);
189+
190+
// Evaluate prolongiation - rhs = B^-1 * N * solCoarse
191+
StdVector prolong(N.rows());
192+
N.multiply(solCoarse, prolong);
193+
B.solve(prolong);
194+
vectorToSolVec(fSim, rhs, prolong);
195+
196+
rhs.beginAssembly();
197+
rhs.endAssembly();
198+
}
199+
200+
//! \brief Returns name of this preconditioner.
201+
const char* getName() const override { return "K-cycle preconditioner"; }
202+
203+
protected:
204+
//! \brief Helper template for copying a equation vector to a DOF vector.
205+
//! \param sim The simulator to use
206+
//! \param V1 The resulting DOF vector
207+
//! \param V2 The initial equation vecto
208+
template <class V1, class V2>
209+
void vectorToSolVec(Sim& sim, V1& dof, const V2& eqn)
210+
{
211+
const int* meqn = sim.getSAM()->getMEQN();
212+
for (size_t i = 1; i <= sim.getPatch(1)->getNoNodes(); ++i) {
213+
int eq = meqn[i-1];
214+
if (eq > 0)
215+
dof(eq) = eqn(i);
216+
}
217+
}
218+
219+
Sim& S1; //!< Preconditioner simulator
220+
Sim& fSim; //!< Solution model
221+
SparseMatrix B, BT; //!< Greville mass matrix for fSim
222+
SparseMatrix N; //!< Basis functions of S1 in greville of fSim
223+
SparseMatrix NT; //!< Stupid no-transpose multiply
224+
};
225+
226+
#endif
227+
228+
#endif

‎Apps/Common/SIMMxV.h‎

Lines changed: 69 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,69 @@
1+
//==============================================================================
2+
//!
3+
//! \file SIMMxV.h
4+
//!
5+
//! \date Nov 28 2012
6+
//!
7+
//! \author Arne Morten Kvarving / SINTEF
8+
//!
9+
//! \brief Use a simulator as a matrix-vector product for iterative solvers.
10+
//!
11+
//==============================================================================
12+
#ifndef SIM_MXV_H_
13+
#define SIM_MXV_H_
14+
15+
#ifdef HAS_PETSC
16+
17+
#include "PETScMatrix.h"
18+
19+
20+
/*!
21+
\brief Template for applying the matrix assembled in a simulator to a vector
22+
for use in iterative solvers.
23+
*/
24+
25+
template<class Sim>
26+
class SIMMxV : public PETScMxV {
27+
public:
28+
//! \brief Default constructor.
29+
//! \param sim The simulator to wrap
30+
SIMMxV(Sim& sim) : S1(sim)
31+
{
32+
S1.getMatrix()->setMxV(this,true);
33+
}
34+
35+
//! \brief Evaluate the matrix-vector product y = A*y
36+
bool evalMxV(Vec& x, Vec& y) override
37+
{
38+
MatMult(S1.getMatrix()->getBlockMatrices()[0], x, y);
39+
return true;
40+
}
41+
42+
//! \brief Set a matrix-free preconditioner.
43+
void setPC(PETScPC* pc)
44+
{
45+
S1.getMatrix()->setPC(pc);
46+
}
47+
48+
//! \brief Solve at time level.
49+
bool solveStep()
50+
{
51+
TimeStep tp;
52+
return S1.solveStep(tp);
53+
}
54+
55+
//! \brief The matrix used for building the preconditioner.
56+
//! \details Defaults to the system matrix
57+
bool evalPC(Mat& P) override
58+
{
59+
MatDuplicate(S1.getMatrix()->getMatrix(), MAT_COPY_VALUES, &P);
60+
return true;
61+
}
62+
63+
protected:
64+
Sim& S1; //!< Reference to wrapped simulator
65+
};
66+
67+
#endif
68+
69+
#endif

0 commit comments

Comments
 (0)