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

Shengyong lua.byhh at gmail.com
Mon Aug 18 00:35:22 CDT 2008


hi, Barry

Thanks for your kindly reply!

Here I report some experiments accoording to your hints.

-pc_type sor -pc_sor_local_symmetric _ksp_type cg failed.
icc(1) also failed.

When using -ksp_monitor_true _residual with icc(0), I have found the precond
residual is osillation around some value about 1.06e-004, the true residual
norm almost staying at 7.128 e-004.  The ||Ae||/||Ax|| also stagnent at
6.3e-005.

After the maximum number of iterations (=2000) reached, the iteration fails.
Even when I set iteration number to 10000, it seems that petsc also fails.

Below is some resiual information near iteration number 2000 when with
ksp_monitor_true_residual.


1990 KSP preconditioned resid norm 1.064720311837e-004 true resid norm
7.1284721
18425e-004 ||Ae||/||Ax|| 6.302380221818e-005

1991 KSP preconditioned resid norm 1.062055494352e-004 true resid norm
7.1281202
17324e-004 ||Ae||/||Ax|| 6.302069101215e-005

1992 KSP preconditioned resid norm 1.061228895583e-004 true resid norm
7.1277740
15661e-004 ||Ae||/||Ax|| 6.301763019565e-005

1993 KSP preconditioned resid norm 1.062165148129e-004 true resid norm
7.1274335
10277e-004 ||Ae||/||Ax|| 6.301461974073e-005

1994 KSP preconditioned resid norm 1.064780917764e-004 true resid norm
7.1270986
90168e-004 ||Ae||/||Ax|| 6.301165955010e-005

1995 KSP preconditioned resid norm 1.068986431546e-004 true resid norm
7.1267695
37853e-004 ||Ae||/||Ax|| 6.300874946922e-005

1996 KSP preconditioned resid norm 1.074687043135e-004 true resid norm
7.1264460
30395e-004 ||Ae||/||Ax|| 6.300588929530e-005

1997 KSP preconditioned resid norm 1.081784773720e-004 true resid norm
7.1261281
40328e-004 ||Ae||/||Ax|| 6.300307878551e-005

1998 KSP preconditioned resid norm 1.090179775109e-004 true resid norm
7.1258158
36472e-004 ||Ae||/||Ax|| 6.300031766417e-005

1999 KSP preconditioned resid norm 1.099771672684e-004 true resid norm
7.1255090
84603e-004 ||Ae||/||Ax|| 6.299760562872e-005

2000 KSP preconditioned resid norm 1.110460758301e-004 true resid norm
7.1252078
48108e-004 ||Ae||/||Ax|| 6.299494235544e-005

On Sun, Aug 17, 2008 at 10:32 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:

>
>   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
>>
>
>


-- 
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/20080818/8357e33b/attachment.htm>


More information about the petsc-users mailing list