strange convergency behavior with petsc (intialy produces good result, suddenly diverges)

[Barry Smith]

    Try using -pc_type sor -pc_sor_local_symmetric -ksp_type cg

    Also try running the original icc one with the additional option - 
ksp_monitor_true_residual, see if funky stuff starts to happen.

    You could also try adding -pc_factor_levels 1 to try ICC(1)  
instead of ICC(0).

    Your code looks fine, I don't see a problem there,

    Barry


    Here is my guess, as the simulation proceeds the variable  
coefficient problem changes enough so that the ICC produces
a badly scaled preconditioner that messes up the iterative method. I  
see this on occasion and don't have a good fix, the shift
positive definite helps sometimes but not always.



On Aug 17, 2008, at 3:36 AM, Shengyong wrote:

> Hi,
>
> I am still struggling to use petsc to solve variable coefficient  
> poisson problems (which orinates from a multi-phase(liquid-gas two  
> phase flow with sharp interface method, the density ratio is 1000,  
> and with surface tension) flow problem) successively.
>
> Initially petsc produces good results with icc pc_type , cg  
> iterative method and with _pc_factor_shift_positive_define flag.  I  
> follow petsc's null space method to constrain the singular property  
> of coefficient matrix. However,  after several time steps good  
> results of simulation, the KSP Residual Norm suddenly reaches to a  
> number greater than 1000. I guess that petsc diverges. I have also  
> tried to swith to other types of methods, e.g. to gmeres, it just  
> behaves similar to cg method.  However, when I use my previous self- 
> coded SOR iterative solver, it produces nice results.  And I have  
> tested the same petsc solver class for a single phase driven cavity  
> flow problem, it also produces nice results.
>
> It seems that I have missed something important setup procedure in  
> my solver class. could anyone to point me the problem ?  I have  
> attached large part of  the code  below:
>
> //in Header
> class CPETScPoissonSolver
> {
>     typedef struct {
>        Vec                  x,b;     //x, b
>        Mat                 A;       //A
>        KSP                ksp;   //Krylov subspace preconditioner
>        PC                   pc;
>        MatNullSpace  nullspace;
>        PetscInt            l, m, n;//L, M, N
>     } UserContext;
>
> public:
>     CPETScPoissonSolver(int argc, char** argv);
>     ~CPETScPoissonSolver(void);
>
>     //........
>
> private:
>
>     //Yale Sparse Matrix format matrix
>      PetscScalar*  A;
>      PetscInt     *   I;
>      PetscInt     *   J;
>
>      // Number of Nonzero Element
>      PetscInt          nnz;
>      //grid step
>      PetscInt          L, M, N;
>      UserContext    userctx;
> private:
>     bool             FirstTime;
> };
>
> //in cpp
> static char helpPetscPoisson[] = "PETSc class Solves a variable  
> Poisson problem with Null Space Method.\n\n";
>
> CPETScPoissonSolver::CPETScPoissonSolver(int argc, char** argv)
> {
>     PetscInitialize(&argc, &argv, (char*)0, helpPetscPoisson);
>     FirstTime=true;
> }
> CPETScPoissonSolver::~CPETScPoissonSolver(void)
> {
>     PetscFinalize();
> }
> ......
> void CPETScPoissonSolver::SetAIJ(PetscScalar *a, PetscInt *i,  
> PetscInt *j, PetscInt Nnz)
> {
>     A= a; I=i; J=j; nnz = Nnz;
> }
>
> PetscErrorCode CPETScPoissonSolver::UserInitializeLinearSolver()
> {
>     PetscErrorCode ierr = 0;
>     PetscInt Num = L*M*N;
>
>     //Since we use successive solvers, so in the second time step we  
> must deallocate the original matrix then setup a new one
>     if(FirstTime==true)
>     {
>         ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,  Num,   
> Num,  I,  J,  A, &userctx.A); CHKERRQ(ierr);
>     }
>     else
>     {
>           FirstTime = false;
>           ierr = MatDestroy(userctx.A);CHKERRQ(ierr);
>           ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,  Num,   
> Num,  I,  J,  A, &userctx.A); CHKERRQ(ierr);
>     }
>
>     if(FirstTime==true)
>     {
>
>         ierr =  
> VecCreateSeqWithArray 
> (PETSC_COMM_SELF,Num,PETSC_NULL,&userctx.b);CHKERRQ(ierr);
>         ierr = VecDuplicate(userctx.b,&userctx.x);CHKERRQ(ierr);
>         ierr = MatNullSpaceCreate(PETSC_COMM_SELF, PETSC_TRUE, 0,  
> PETSC_NULL, &userctx.nullspace); CHKERRQ(ierr);
>         ierr = KSPCreate(PETSC_COMM_SELF,&userctx.ksp);CHKERRQ(ierr);
>         /*Set Null Space for KSP*/
>         ierr = KSPSetNullSpace(userctx.ksp,  
> userctx.nullspace);CHKERRQ(ierr);
>     }
>     return 0;
> }
>
>
> PetscErrorCode CPETScPoissonSolver::UserSetBX(PetscScalar *x,  
> PetscScalar *b)
> {
>       PetscErrorCode ierr ;
>       //below code we must set it every time step
>       ierr = VecPlaceArray(userctx.x,x);CHKERRQ(ierr);
>       ierr = VecPlaceArray(userctx.b,b);CHKERRQ(ierr);
>       ierr =  MatNullSpaceRemove(userctx.nullspace,userctx.b,  
> PETSC_NULL);CHKERRQ(ierr);
>       return 0;
> }
>
> PetscInt CPETScPoissonSolver::UserSolve()
> {
>       PetscErrorCode ierr;
>       ierr =  
> KSPSetOperators 
> (userctx.ksp,userctx.A,userctx.A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
>       ierr = KSPSetType(userctx.ksp, KSPCG);
>       ierr = KSPSetInitialGuessNonzero(userctx.ksp, PETSC_TRUE);
>       ierr = KSPGetPC(userctx.ksp,&userctx.pc);CHKERRQ(ierr);
>       ierr = PCSetType(userctx.pc,PCICC);CHKERRQ(ierr);
>       ierr = PCFactorSetShiftPd(userctx.pc, PETSC_TRUE);
>       ierr = KSPSetTolerances(userctx.ksp, 
> 1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,2000);
>       ierr = KSPSetFromOptions(userctx.ksp);CHKERRQ(ierr);
>       /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  
> - - - -
>                           Solve the linear system
>          - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  
> - - - - */
>       ierr = KSPSolve(userctx.ksp,userctx.b,userctx.x);CHKERRQ(ierr);
>
>       ierr = VecResetArray(userctx.x);CHKERRQ(ierr);
>       ierr = VecResetArray(userctx.b);CHKERRQ(ierr);
>
>       return 0;
> }
>
> PetscErrorCode CPETScPoissonSolver::ReleaseMem()
> {
>       PetscErrorCode ierr;
>       ierr = KSPDestroy(userctx.ksp);CHKERRQ(ierr);
>       ierr = VecDestroy(userctx.x);     CHKERRQ(ierr);
>       ierr = VecDestroy(userctx.b);     CHKERRQ(ierr);
>       ierr = MatDestroy(userctx.A);    CHKERRQ(ierr);
>       ierr = MatNullSpaceDestroy(userctx.nullspace); CHKERRQ(ierr);
>       return 0;
> }
>
> Thanks very much!
>
>
> -- 
> Pang Shengyong
> Solidification Simulation Lab,
> State Key Lab of Mould & Die Technology,
> Huazhong Univ. of Sci. & Tech. China