Thanks Barry. Yes that was my problem now if I run with the same down and up number of iterations I see the in -ksp_view output that smoother up is the same as smoother down.
I think I figure it out how to set up different smoothers up and down but use the same ASM Preconditioner, which is more or less what Lorenzo suggested. Thanks again Eugenio ________________________________________ From: Barry Smith [[email protected]] Sent: Thursday, August 20, 2015 11:54 PM To: Aulisa, Eugenio Cc: PETSc list Subject: Re: [petsc-users] GMRES -> PCMG -> PCASM pre- post- smoother Ahhh, void AsmPetscLinearEquationSolver::MGsolve ( const bool ksp_clean , const unsigned &npre, const unsigned &npost ) { if ( ksp_clean ) { PetscMatrix* KKp = static_cast< PetscMatrix* > ( _KK ); Mat KK = KKp->mat(); KSPSetOperators ( _ksp, KK, _Pmat ); KSPSetTolerances ( _ksp, _rtol, _abstol, _dtol, _maxits ); KSPSetFromOptions ( _ksp ); PC pcMG; KSPGetPC(_ksp, &pcMG); PCMGSetNumberSmoothDown(pcMG, npre); PCMGSetNumberSmoothUp(pcMG, npost); } PetscErrorCode PCMGSetNumberSmoothDown(PC pc,PetscInt n) { PC_MG *mg = (PC_MG*)pc->data; PC_MG_Levels **mglevels = mg->levels; PetscErrorCode ierr; PetscInt i,levels; PetscFunctionBegin; PetscValidHeaderSpecific(pc,PC_CLASSID,1); if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling"); PetscValidLogicalCollectiveInt(pc,n,2); levels = mglevels[0]->levels; for (i=1; i<levels; i++) { /* make sure smoother up and down are different */ ierr = PCMGGetSmootherUp(pc,i,NULL);CHKERRQ(ierr); ierr = KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);CHKERRQ(ierr); mg->default_smoothd = n; } PetscFunctionReturn(0); } PetscErrorCode PCMGGetSmootherUp(PC pc,PetscInt l,KSP *ksp) { PC_MG *mg = (PC_MG*)pc->data; PC_MG_Levels **mglevels = mg->levels; PetscErrorCode ierr; const char *prefix; MPI_Comm comm; PetscFunctionBegin; PetscValidHeaderSpecific(pc,PC_CLASSID,1); /* This is called only if user wants a different pre-smoother from post. Thus we check if a different one has already been allocated, if not we allocate it. */ if (!l) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"There is no such thing as a up smoother on the coarse grid"); if (mglevels[l]->smoothu == mglevels[l]->smoothd) { KSPType ksptype; PCType pctype; PC ipc; PetscReal rtol,abstol,dtol; PetscInt maxits; KSPNormType normtype; ierr = PetscObjectGetComm((PetscObject)mglevels[l]->smoothd,&comm);CHKERRQ(ierr); ierr = KSPGetOptionsPrefix(mglevels[l]->smoothd,&prefix);CHKERRQ(ierr); ierr = KSPGetTolerances(mglevels[l]->smoothd,&rtol,&abstol,&dtol,&maxits);CHKERRQ(ierr); ierr = KSPGetType(mglevels[l]->smoothd,&ksptype);CHKERRQ(ierr); ierr = KSPGetNormType(mglevels[l]->smoothd,&normtype);CHKERRQ(ierr); ierr = KSPGetPC(mglevels[l]->smoothd,&ipc);CHKERRQ(ierr); ierr = PCGetType(ipc,&pctype);CHKERRQ(ierr); ierr = KSPCreate(comm,&mglevels[l]->smoothu);CHKERRQ(ierr); ierr = KSPSetErrorIfNotConverged(mglevels[l]->smoothu,pc->erroriffailure);CHKERRQ(ierr); ierr = PetscObjectIncrementTabLevel((PetscObject)mglevels[l]->smoothu,(PetscObject)pc,mglevels[0]->levels-l);CHKERRQ(ierr); ierr = KSPSetOptionsPrefix(mglevels[l]->smoothu,prefix);CHKERRQ(ierr); ierr = KSPSetTolerances(mglevels[l]->smoothu,rtol,abstol,dtol,maxits);CHKERRQ(ierr); ierr = KSPSetType(mglevels[l]->smoothu,ksptype);CHKERRQ(ierr); ierr = KSPSetNormType(mglevels[l]->smoothu,normtype);CHKERRQ(ierr); ierr = KSPSetConvergenceTest(mglevels[l]->smoothu,KSPConvergedSkip,NULL,NULL);CHKERRQ(ierr); ierr = KSPGetPC(mglevels[l]->smoothu,&ipc);CHKERRQ(ierr); ierr = PCSetType(ipc,pctype);CHKERRQ(ierr); ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)mglevels[l]->smoothu);CHKERRQ(ierr); } if (ksp) *ksp = mglevels[l]->smoothu; PetscFunctionReturn(0); } As soon as you set both the up and down number of iterations it causes a duplication of the current smoother with some options preserved but others not (we don't have a KSPDuplicate() that duplicates everything). So if you are fine with the number of pre and post smooths the same just don't set both PCMGSetNumberSmoothDown(pcMG, npre); PCMGSetNumberSmoothUp(pcMG, npost); if you want them to be different you can share the same PC between the two (which has the overlapping matrices in it) but you cannot share the same KSP. I can tell you how to do that but suggest it is simpler just to have the same number of pre and post smooths Barry > On Aug 20, 2015, at 6:51 AM, Aulisa, Eugenio <[email protected]> wrote: > > Hi Barry, > > Thanks for your answer. > > I run my applications with no command line, and I do not think I changed any > PETSC_OPTIONS, > at least not voluntarily. > > For the source it is available on > https://github.com/NumPDEClassTTU/femus > but it is part of a much larger library and > I do not think any of you want to install and run it > just to find what I messed up. > > In any case, if you just want to look at the source code > where I set up the level smoother it is in > > https://github.com/NumPDEClassTTU/femus/blob/master/src/algebra/AsmPetscLinearEquationSolver.cpp > > line 400 > > void AsmPetscLinearEquationSolver::MGsetLevels ( > LinearEquationSolver *LinSolver, const unsigned &level, const unsigned > &levelMax, > const vector <unsigned> &variable_to_be_solved, SparseMatrix* PP, > SparseMatrix* RR ){ > > Be aware, that even if it seams that this takes care of the coarse level it > is not. > The coarse level smoother is set some where else. > > Thanks, > Eugenio > > ________________________________________ > From: Barry Smith [[email protected]] > Sent: Thursday, August 20, 2015 2:37 AM > To: Aulisa, Eugenio > Cc: [email protected] > Subject: Re: [petsc-users] GMRES -> PCMG -> PCASM pre- post- smoother > > What you describe is not the expected behavior. I expected exactly the > result that you expected. > > Do you perhaps have some PETSc options around that may be changing the > post-smoother? On the command line or in the file petscrc or in the > environmental variable PETSC_OPTIONS? Can you send us some code that we > could run that reproduces the problem? > > Barry > >> On Aug 19, 2015, at 9:26 PM, Aulisa, Eugenio <[email protected]> wrote: >> >> Hi, >> >> I am solving an iteration of >> >> GMRES -> PCMG -> PCASM >> >> where I build my particular ASM domain decomposition. >> >> In setting the PCMG I would like at each level >> to use the same pre- and post-smoother >> and for this reason I am using >> ... >> PCMGGetSmoother ( pcMG, level , &subksp ); >> >> to extract and set at each level the ksp object. >> >> In setting PCASM then I use >> ... >> KSPGetPC ( subksp, &subpc ); >> PCSetType ( subpc, PCASM ); >> ... >> and then set my own decomposition >> ... >> PCASMSetLocalSubdomains(subpc,_is_loc_idx.size(),&_is_ovl[0],&_is_loc[0]); >> ... >> >> Now everything compiles, and runs with no memory leakage, >> but I do not get the expected convergence. >> >> When I checked the output of -ksp_view, I saw something that puzzled me: >> at each level >0, while in the MG pre-smoother the ASM domain decomposition >> is the one that I set, for example with 4 processes I get >> >>>>>>>>>>>>>>>>>>>>> >> ... >> Down solver (pre-smoother) on level 2 ------------------------------- >> KSP Object: (level-2) 4 MPI processes >> type: gmres >> GMRES: restart=30, using Classical (unmodified) Gram-Schmidt >> Orthogonalization with no iterative refinement >> GMRES: happy breakdown tolerance 1e-30 >> maximum iterations=1 >> using preconditioner applied to right hand side for initial guess >> tolerances: relative=1e-12, absolute=1e-20, divergence=1e+50 >> left preconditioning >> using nonzero initial guess >> using NONE norm type for convergence test >> PC Object: (level-2) 4 MPI processes >> type: asm >> Additive Schwarz: total subdomain blocks = 198, amount of overlap = 0 >> Additive Schwarz: restriction/interpolation type - RESTRICT >> [0] number of local blocks = 52 >> [1] number of local blocks = 48 >> [2] number of local blocks = 48 >> [3] number of local blocks = 50 >> Local solve info for each block is in the following KSP and PC objects: >> - - - - - - - - - - - - - - - - - - >> ... >>>>>>>>>>>>> >> >> >> in the post-smoother I have the default ASM decomposition with overlapping >> 1: >> >> >>>>>>>>>>>>> >> ... >> Up solver (post-smoother) on level 2 ------------------------------- >> KSP Object: (level-2) 4 MPI processes >> type: gmres >> GMRES: restart=30, using Classical (unmodified) Gram-Schmidt >> Orthogonalization with no iterative refinement >> GMRES: happy breakdown tolerance 1e-30 >> maximum iterations=2 >> tolerances: relative=1e-12, absolute=1e-20, divergence=1e+50 >> left preconditioning >> using nonzero initial guess >> using NONE norm type for convergence test >> PC Object: (level-2) 4 MPI processes >> type: asm >> Additive Schwarz: total subdomain blocks = 4, amount of overlap = 1 >> Additive Schwarz: restriction/interpolation type - RESTRICT >> Local solve is same for all blocks, in the following KSP and PC >> objects: >> KSP Object: (level-2sub_) 1 MPI processes >> type: preonly >> maximum iterations=10000, initial guess is zero >> tolerances: relative=1e-05, absolute=1e-50, divergence=10000 >> left preconditioning >> ... >>>>>>>>>>>>>>> >> %%%%%%%%%%%%%%%%%%%%%%%% >> >> So it seams that by using >> >> PCMGGetSmoother ( pcMG, level , &subksp ); >> >> I was capable to set both the pre- and post- smoothers to be PCASM >> but everything I did after that applied only to the >> pre-smoother, while the post-smoother got the default PCASM options. >> >> I know that I can use >> PCMGGetSmootherDown and PCMGGetSmootherUp, but that would >> probably double the memory allocation and the computational time in the ASM. >> >> Is there any way I can just use PCMGGetSmoother >> and use the same PCASM in the pre- and post- smoother? >> >> I hope I was clear enough. >> >> Thanks a lot for your help, >> Eugenio >> >> >> >
