Skip to content

Commit

Permalink
Moved _PP and _RR form LinearEquation to ImplicitLinearSystem
Browse files Browse the repository at this point in the history
  • Loading branch information
@eaulisa
eaulisa committed Mar 20, 2015
1 parent 138f923 commit 7388fe1
Show file tree
Hide file tree
Showing 6 changed files with 39 additions and 30 deletions.
1 change: 0 additions & 1 deletion src/algebra/LinearEquation.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,6 @@ class LinearEquation : public ParallelObject {
// member data
Mesh *_msh;
NumericVector *_EPS, *_EPSC, *_RES, *_RESC;
//SparseMatrix *_KK, *_CC;
SparseMatrix *_KK, *_PP,*_RR, *_CC;
vector < vector <unsigned> > KKoffset;
vector < unsigned > KKghostsize;
Expand Down
39 changes: 24 additions & 15 deletions src/equations/LinearImplicitSystem.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,10 @@ void LinearImplicitSystem::clear() {
for (unsigned ig=0; ig<_gridn; ig++) {
_LinSolver[ig]->DeletePde();
delete _LinSolver[ig];
}
if(_PP[ig]) delete _PP[ig];
if(_RR[ig]) delete _RR[ig];
}

_NSchurVar_test=0;
_numblock_test=0;
_numblock_all_test=0;
Expand Down Expand Up @@ -87,7 +90,13 @@ void LinearImplicitSystem::init() {
_LinSolver[i]->InitPde(_SolSystemPdeIndex,_ml_sol->GetSolType(),
_ml_sol->GetSolName(),&_solution[i]->_Bdc,_gridr,_gridn,_SparsityPattern);
}


_PP.resize(_gridn);
_RR.resize(_gridn);
for(unsigned i=0;i<_gridn;i++){
_PP[i]=NULL;
_RR[i]=NULL;
}
for (unsigned ig=1; ig<_gridn; ig++) {
BuildProlongatorMatrix(ig);
}
Expand Down Expand Up @@ -394,32 +403,32 @@ void LinearImplicitSystem::Restrictor(const unsigned &gridf, const unsigned &gri
if (gridf>=_gridr) { //_gridr
if (!_LinSolver[gridf-1]->_CC_flag) {
_LinSolver[gridf-1]->_CC_flag=1;
_LinSolver[gridf-1]->_CC->matrix_PtAP(*_LinSolver[gridf]->_PP,*_LinSolver[gridf]->_KK,!matrix_reuse);
_LinSolver[gridf-1]->_CC->matrix_PtAP(*_PP[gridf],*_LinSolver[gridf]->_KK,!matrix_reuse);
}
else{
_LinSolver[gridf-1]->_CC->matrix_PtAP(*_LinSolver[gridf]->_PP,*_LinSolver[gridf]->_KK,matrix_reuse);
_LinSolver[gridf-1]->_CC->matrix_PtAP(*_PP[gridf],*_LinSolver[gridf]->_KK,matrix_reuse);
}
_LinSolver[gridf-1u]->_KK->matrix_add(1.,*_LinSolver[gridf-1u]->_CC,"different_nonzero_pattern");
}
else { //Projection of the Matrix on the lower level
if (non_linear_iteration==0 && ( full_cycle*(gridf==gridn-1u) || !full_cycle )) {
_LinSolver[gridf-1]->_KK->matrix_PtAP(*_LinSolver[gridf]->_PP,*_LinSolver[gridf]->_KK,!matrix_reuse);
_LinSolver[gridf-1]->_KK->matrix_PtAP(*_PP[gridf],*_LinSolver[gridf]->_KK,!matrix_reuse);
}
else{
_LinSolver[gridf-1]->_KK->matrix_PtAP(*_LinSolver[gridf]->_PP,*_LinSolver[gridf]->_KK,matrix_reuse);
_LinSolver[gridf-1]->_KK->matrix_PtAP(*_PP[gridf],*_LinSolver[gridf]->_KK,matrix_reuse);
}
}
}

_LinSolver[gridf-1u]->_RESC->matrix_mult_transpose(*_LinSolver[gridf]->_RES, *_LinSolver[gridf]->_PP);
_LinSolver[gridf-1u]->_RESC->matrix_mult_transpose(*_LinSolver[gridf]->_RES, *_PP[gridf]);
*_LinSolver[gridf-1u]->_RES += *_LinSolver[gridf-1u]->_RESC;

}

// *******************************************************
void LinearImplicitSystem::Prolongator(const unsigned &gridf) {

_LinSolver[gridf]->_EPSC->matrix_mult(*_LinSolver[gridf-1]->_EPS,*_LinSolver[gridf]->_PP);
_LinSolver[gridf]->_EPSC->matrix_mult(*_LinSolver[gridf-1]->_EPS,*_PP[gridf]);
_LinSolver[gridf]->UpdateResidual();
_LinSolver[gridf]->SumEpsCToEps();
}
Expand Down Expand Up @@ -498,8 +507,8 @@ void LinearImplicitSystem::BuildProlongatorMatrix(unsigned gridf) {
delete NNZ_d;
delete NNZ_o;

LinSolf->_PP = SparseMatrix::build().release();
LinSolf->_PP->init(nf,nc,nf_loc,nc_loc,nnz_d,nnz_o);
_PP[gridf] = SparseMatrix::build().release();
_PP[gridf]->init(nf,nc,nf_loc,nc_loc,nnz_d,nnz_o);

for (unsigned k=0; k<_SolSystemPdeIndex.size(); k++) {
unsigned SolIndex = _SolSystemPdeIndex[k];
Expand All @@ -510,12 +519,12 @@ void LinearImplicitSystem::BuildProlongatorMatrix(unsigned gridf) {
unsigned iel = mshc->IS_Mts2Gmt_elem[iel_mts];
if(mshc->el->GetRefinedElementIndex(iel)){ //only if the coarse element has been refined
short unsigned ielt=mshc->el->GetElementType(iel);
mshc->_finiteElement[ielt][SolType]->BuildProlongation(*LinSolf,*LinSolc,iel,LinSolf->_PP,SolIndex,k);
mshc->_finiteElement[ielt][SolType]->BuildProlongation(*LinSolf,*LinSolc,iel,_PP[gridf],SolIndex,k);
}
}
}
}
LinSolf->_PP->close();
_PP[gridf]->close();
}


Expand Down Expand Up @@ -683,7 +692,7 @@ void LinearImplicitSystem::MGSolve(double Eps1, // tolerance for the li

for (uint Level = 1; Level < GetGridn(); Level++) {

_LinSolver[Level]->_EPS->matrix_mult(*_LinSolver[Level-1]->_EPS,*_LinSolver[Level]->_PP); //**** project the solution
_LinSolver[Level]->_EPS->matrix_mult(*_LinSolver[Level-1]->_EPS,*_PP[Level]); //**** project the solution

res_fine = MGStep(Level,Eps1,MaxIter,Gamma,Nc_pre,Nc_coarse,Nc_post);

Expand Down Expand Up @@ -808,7 +817,7 @@ double LinearImplicitSystem::MGStep(int Level, // Level

/// std::cout << ">>>>>>>> BEGIN ONE DESCENT >>>>>>>>>>"<< std::endl;

_LinSolver[Level-1]->_RESC->matrix_mult(*_LinSolver[Level]->_RES,*_LinSolver[Level-1]->_RR);//****** restrict the residual from the finer grid ( new rhs )
_LinSolver[Level-1]->_RESC->matrix_mult(*_LinSolver[Level]->_RES,*_RR[Level-1]);//****** restrict the residual from the finer grid ( new rhs )
_LinSolver[Level-1]->_EPS->close(); //initial value of x for the presmoothing iterations
_LinSolver[Level-1]->_EPS->zero(); //initial value of x for the presmoothing iterations

Expand All @@ -823,7 +832,7 @@ double LinearImplicitSystem::MGStep(int Level, // Level
std::cout << "BEFORE PROL Level " << Level << " res linfty " << _LinSolver[Level-1]->_RES->linfty_norm() << " res l2 " << _LinSolver[Level-1]->_RES->l2_norm() << std::endl;
#endif

_LinSolver[Level]->_RES->matrix_mult(*_LinSolver[Level-1]->_EPS,*_LinSolver[Level]->_PP);//******** project the dx from the coarser grid
_LinSolver[Level]->_RES->matrix_mult(*_LinSolver[Level-1]->_EPS,*_PP[Level]);//******** project the dx from the coarser grid
#ifdef DEFAULT_PRINT_CONV
_LinSolver[Level]->_RES->close();
//here, the _res contains the prolongation of dx, so the new x is x_old + P dx
Expand Down
3 changes: 2 additions & 1 deletion src/equations/LinearImplicitSystem.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ class LinearImplicitSystem : public ImplicitSystem {
* like PETSc or LASPACK. Up to now also for the nonlinear case we use linear_solvers, in future we will add the nonlinear solver
*/
vector < LinearEquationSolver*> _LinSolver;
vector < SparseMatrix* > _PP, _RR;

/** Set the max number of linear iterationsfor solving Ax=b */
void SetMaxNumberOfLinearIterations(unsigned int max_lin_it) {
_n_max_linear_iterations = max_lin_it;
Expand Down Expand Up @@ -170,6 +170,7 @@ class LinearImplicitSystem : public ImplicitSystem {

protected:

vector < SparseMatrix* > _PP, _RR;
/** Create the Prolongator matrix for the Multigrid solver */
void Prolongator(const unsigned &gridf);

Expand Down
22 changes: 11 additions & 11 deletions src/equations/MonolithicFSINonLinearImplicitSystem.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,24 +65,24 @@ void MonolithicFSINonLinearImplicitSystem::Restrictor(const unsigned &gridf, con
if (gridf>=_gridr) { //_gridr
if (!LinSolc->_CC_flag) {
LinSolc->_CC_flag=1;
LinSolc->_CC->matrix_ABC(*LinSolf->_RR,*LinSolf->_KK,*LinSolf->_PP,!matrix_reuse);
LinSolc->_CC->matrix_ABC(*_RR[gridf],*LinSolf->_KK,*_PP[gridf],!matrix_reuse);
}
else{
LinSolc->_CC->matrix_ABC(*LinSolf->_RR,*LinSolf->_KK,*LinSolf->_PP,matrix_reuse);
LinSolc->_CC->matrix_ABC(*_RR[gridf],*LinSolf->_KK,*_PP[gridf],matrix_reuse);
}
LinSolc->_KK->matrix_add(1.,*LinSolc->_CC,"different_nonzero_pattern");
}
else { //Projection of the Matrix on the lower level
if (non_linear_iteration==0 && ( full_cycle*(gridf==gridn-1u) || !full_cycle )) {
LinSolc->_KK->matrix_ABC(*LinSolf->_RR,*LinSolf->_KK,*LinSolf->_PP,!matrix_reuse);
LinSolc->_KK->matrix_ABC(*_RR[gridf],*LinSolf->_KK,*_PP[gridf],!matrix_reuse);
}
else{
LinSolc->_KK->matrix_ABC(*LinSolf->_RR,*LinSolf->_KK,*LinSolf->_PP,matrix_reuse);
LinSolc->_KK->matrix_ABC(*_RR[gridf],*LinSolf->_KK,*_PP[gridf],matrix_reuse);
}
}
}

LinSolc->_RESC->matrix_mult(*LinSolf->_RES, *LinSolf->_RR);
LinSolc->_RESC->matrix_mult(*LinSolf->_RES, *_RR[gridf]);
*LinSolc->_RES += *LinSolc->_RESC;

}
Expand Down Expand Up @@ -145,8 +145,8 @@ void MonolithicFSINonLinearImplicitSystem::BuildProlongatorMatrix(unsigned gridf
delete NNZ_d;
delete NNZ_o;

LinSolf->_PP = SparseMatrix::build().release();
LinSolf->_PP->init(nf,nc,nf_loc,nc_loc,nnz_d, nnz_o);
_PP[gridf] = SparseMatrix::build().release();
_PP[gridf]->init(nf,nc,nf_loc,nc_loc,nnz_d, nnz_o);

SparseMatrix *RRt;
RRt = SparseMatrix::build().release();
Expand Down Expand Up @@ -179,17 +179,17 @@ void MonolithicFSINonLinearImplicitSystem::BuildProlongatorMatrix(unsigned gridf
else{
mshc->_finiteElement[ielt][SolType]->BuildProlongation(*LinSolf,*LinSolc,iel, RRt,SolIndex,k);
}
mshc->_finiteElement[ielt][SolType]->BuildProlongation(*LinSolf,*LinSolc,iel, LinSolf->_PP,SolIndex,k);
mshc->_finiteElement[ielt][SolType]->BuildProlongation(*LinSolf,*LinSolc,iel, _PP[gridf],SolIndex,k);
}
}
}
}

LinSolf->_PP->close();
_PP[gridf]->close();
RRt->close();

LinSolf->_RR = SparseMatrix::build().release();
RRt->get_transpose( *LinSolf->_RR);
_RR[gridf] = SparseMatrix::build().release();
RRt->get_transpose( *_RR[gridf]);
delete RRt;

}
Expand Down
2 changes: 1 addition & 1 deletion src/solution/GMVWriter.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -402,7 +402,7 @@ void GMVWriter::Pwrite(const std::string output_path, const char order[], const

if(n_processors()==1) { // IF SERIAL
Mysol[ig]->init(_ml_mesh->GetLevel(ig)->MetisOffset[index][_nprocs],
_ml_mesh->GetLevel(ig)->own_size[index][_nprocs],false,SERIAL);
_ml_mesh->GetLevel(ig)->MetisOffset[index][_nprocs],false,SERIAL);
}
else{ // IF PARALLEL
Mysol[ig]->init(_ml_mesh->GetLevel(ig)->MetisOffset[index][_nprocs],
Expand Down
2 changes: 1 addition & 1 deletion src/solution/VTKWriter.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -174,7 +174,7 @@ void VTKWriter::Pwrite(const std::string output_path, const char order[], const

if(n_processors()==1) { // IF SERIAL
mysol[ig]->init(_ml_mesh->GetLevel(ig)->MetisOffset[index][_nprocs],
_ml_mesh->GetLevel(ig)->own_size[index][_nprocs],false,SERIAL);
_ml_mesh->GetLevel(ig)->MetisOffset[index][_nprocs],false,SERIAL);
}
else{ // IF PARALLEL
mysol[ig]->init(_ml_mesh->GetLevel(ig)->MetisOffset[index][_nprocs],
Expand Down

0 comments on commit 7388fe1

Please sign in to comment.