strange convergency behavior with petsc (intialy produces good result, suddenly diverges)
Shengyong
lua.byhh at gmail.com
Tue Aug 19 01:07:42 CDT 2008
If do not use PETSc or laspack,the result is fine. I have comparied with
some literature results.
Maybe there are some bugs in my coeffient matrix generation code. I will try
to fix it.
Thanks.
On Tue, Aug 19, 2008 at 3:23 AM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
> 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
>>
>
>
--
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/d956ea72/attachment.htm>
More information about the petsc-users
mailing list