strange convergency behavior with petsc (intialy produces good result, suddenly diverges)
Shengyong
lua.byhh at gmail.com
Mon Aug 18 12:02:18 CDT 2008
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20080819/dad2daf0/attachment.htm>
More information about the petsc-users
mailing list