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

[Shengyong]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20080817/5c312e40/attachment.htm>