strange convergency behavior with petsc (intialy produces good result, suddenly diverges)
Barry Smith
bsmith at mcs.anl.gov
Mon Aug 18 14:23:38 CDT 2008
Again, what happens with -pc_type lu? Does the simulation run for
any number of time steps with correct answers?
Barry
On Aug 18, 2008, at 12:02 PM, Shengyong wrote:
> Hi, Barry
>
> Thanks for your help!
>
> What I mean is that both of your suggested method converge for a
> while then stagnant. Once stagnation happens, then after maximum
> iteration number reaches, PETSc fill the solution with exactly
> zero. So error occurs in the next time step. I am a newbie and can
> not explain why this happens since I have set the initial guess is
> not zero in my code. E.g, I have called this function
> KSPSetInitialGuessNonzero(userctx.ksp, PETSC_TRUE);
>
> For the stagnation, it seems that the rounding error dominates
> during iteration in some special situations since the coefficient
> Matrix A is very ill-conditioned(the Cond number over 10^6 in some
> cases). Although with null space method by project B vector to R(A)
> will remede the stagnant phenomena. Actually it does not. I have
> read a paper in J. Comp. Appl. Math. in year of 1988 :
> 'preconditioned conjugate gradient for solving singular systems'.
> The author suggests not only project B to R(A) space, but also X(i
> +1) (solution vector during iteration) and r(i+1) ( residual vector
> during iteration) to R(A*) (A* refers to preconditioned A). I am not
> clear about how icc and cg method are implemented in PETSc. Above is
> just my not so mature guess.
>
> My simulation code could run unlimited time steps with self-coded
> SOR method. The SOR method is implemented as matrix free method. In
> this method, the Neumann BC is not directly discretized in
> equations, but this way is adopted: assuming P[0] is BC, then P[0]
> cell does not particpate in any iteration, but rather after one step
> iteration ends, I set P[0] = P[1] while P[1] cell is inner fluid
> cell. While with PETSc method, I directly discrete the Neumann BC
> into boundary cells' equation. I think this difference does not
> matter the converge stuff at all since the initially severay time
> steps PETSc solver gives nice results.
>
> Best Regards,
>
>
>
> On Mon, Aug 18, 2008 at 11:24 PM, Barry Smith <bsmith at mcs.anl.gov>
> wrote:
>
> On Aug 18, 2008, at 12:35 AM, Shengyong wrote:
>
> 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.
>
> What does failed mean? Converged for a while then stagnated?
> Diverged (the residual norm got larger and larger)? Never converged
> at all?
>
> icc(1) also failed.
>
> What does failed mean?
>
> I suggest running with -pc_type lu (keep the -ksp_type gmres) and
> see what happens? Does you simulation run fine for an unlimited number
> of timesteps? If not, how does it fail?
>
> Barry
>
>
>
>
> 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
>
>
>
>
> --
> Pang Shengyong
> Solidification Simulation Lab,
> State Key Lab of Mould & Die Technology,
> Huazhong Univ. of Sci. & Tech. China
More information about the petsc-users
mailing list