<br><br><div class="gmail_quote">On Mon, Aug 18, 2008 at 11:24 PM, Barry Smith <span dir="ltr">&lt;<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>&gt;</span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<div class="Ih2E3d"><br>On Aug 18, 2008, at 12:35 AM, Shengyong wrote:<br><br>
<blockquote class="gmail_quote" style="PADDING-LEFT: 1ex; MARGIN: 0px 0px 0px 0.8ex; BORDER-LEFT: #ccc 1px solid">hi, Barry<br><br>Thanks for your kindly reply!<br><br>Here I report some experiments accoording to your hints.<br>
<br>-pc_type sor -pc_sor_local_symmetric _ksp_type cg failed.<br></blockquote><br></div>&nbsp; What does failed mean? Converged for a while then stagnated? Diverged (the residual norm got larger and larger)? Never converged at all?<br>

<blockquote class="gmail_quote" style="PADDING-LEFT: 1ex; MARGIN: 0px 0px 0px 0.8ex; BORDER-LEFT: #ccc 1px solid"><br>icc(1) also failed.<br></blockquote><br>&nbsp; What does failed mean?<br><br>&nbsp; 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<br>
of timesteps? If not, how does it fail?<br><font color="#888888"><br>&nbsp; Barry</font>
<div>
<div></div>
<div class="Wj3C7c"><br><br>
<blockquote class="gmail_quote" style="PADDING-LEFT: 1ex; MARGIN: 0px 0px 0px 0.8ex; BORDER-LEFT: #ccc 1px solid"><br><br>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. &nbsp;The ||Ae||/||Ax|| also stagnent at 6.3e-005.<br>
<br>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.<br><br>Below is some resiual information near iteration number 2000 when with ksp_monitor_true_residual.<br>
<br><br>1990 KSP preconditioned resid norm 1.064720311837e-004 true resid norm 7.1284721<br>18425e-004 ||Ae||/||Ax|| 6.302380221818e-005<br><br>1991 KSP preconditioned resid norm 1.062055494352e-004 true resid norm 7.1281202<br>
17324e-004 ||Ae||/||Ax|| 6.302069101215e-005<br><br>1992 KSP preconditioned resid norm 1.061228895583e-004 true resid norm 7.1277740<br>15661e-004 ||Ae||/||Ax|| 6.301763019565e-005<br><br>1993 KSP preconditioned resid norm 1.062165148129e-004 true resid norm 7.1274335<br>
10277e-004 ||Ae||/||Ax|| 6.301461974073e-005<br><br>1994 KSP preconditioned resid norm 1.064780917764e-004 true resid norm 7.1270986<br>90168e-004 ||Ae||/||Ax|| 6.301165955010e-005<br><br>1995 KSP preconditioned resid norm 1.068986431546e-004 true resid norm 7.1267695<br>
37853e-004 ||Ae||/||Ax|| 6.300874946922e-005<br><br>1996 KSP preconditioned resid norm 1.074687043135e-004 true resid norm 7.1264460<br>30395e-004 ||Ae||/||Ax|| 6.300588929530e-005<br><br>1997 KSP preconditioned resid norm 1.081784773720e-004 true resid norm 7.1261281<br>
40328e-004 ||Ae||/||Ax|| 6.300307878551e-005<br><br>1998 KSP preconditioned resid norm 1.090179775109e-004 true resid norm 7.1258158<br>36472e-004 ||Ae||/||Ax|| 6.300031766417e-005<br><br>1999 KSP preconditioned resid norm 1.099771672684e-004 true resid norm 7.1255090<br>
84603e-004 ||Ae||/||Ax|| 6.299760562872e-005<br><br>2000 KSP preconditioned resid norm 1.110460758301e-004 true resid norm 7.1252078<br>48108e-004 ||Ae||/||Ax|| 6.299494235544e-005<br><br>On Sun, Aug 17, 2008 at 10:32 PM, Barry Smith &lt;<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>&gt; wrote:<br>
<br>&nbsp;Try using -pc_type sor -pc_sor_local_symmetric -ksp_type cg<br><br>&nbsp;Also try running the original icc one with the additional option -ksp_monitor_true_residual, see if funky stuff starts to happen.<br><br>&nbsp;You could also try adding -pc_factor_levels 1 to try ICC(1) instead of ICC(0).<br>
<br>&nbsp;Your code looks fine, I don&#39;t see a problem there,<br><br>&nbsp;Barry<br><br><br>&nbsp;Here is my guess, as the simulation proceeds the variable coefficient problem changes enough so that the ICC produces<br>a badly scaled preconditioner that messes up the iterative method. I see this on occasion and don&#39;t have a good fix, the shift<br>
positive definite helps sometimes but not always.<br><br><br><br><br>On Aug 17, 2008, at 3:36 AM, Shengyong wrote:<br><br>Hi,<br><br>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.<br>
<br>Initially petsc produces good results with icc pc_type , cg iterative method and with _pc_factor_shift_positive_define flag. &nbsp;I follow petsc&#39;s null space method to constrain the singular property of coefficient matrix. However, &nbsp;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. &nbsp;However, when I use my previous self-coded SOR iterative solver, it produces nice results. &nbsp;And I have tested the same petsc solver class for a single phase driven cavity flow problem, it also produces nice results.<br>
<br>It seems that I have missed something important setup procedure in my solver class. could anyone to point me the problem ? &nbsp;I have attached large part of &nbsp;the code &nbsp;below:<br><br>//in Header<br>class CPETScPoissonSolver<br>
{<br>&nbsp; typedef struct {<br>&nbsp; &nbsp; &nbsp;Vec &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;x,b; &nbsp; &nbsp; //x, b<br>&nbsp; &nbsp; &nbsp;Mat &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; A; &nbsp; &nbsp; &nbsp; //A<br>&nbsp; &nbsp; &nbsp;KSP &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;ksp; &nbsp; //Krylov subspace preconditioner<br>&nbsp; &nbsp; &nbsp;PC &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; pc;<br>&nbsp; &nbsp; &nbsp;MatNullSpace &nbsp;nullspace;<br>
&nbsp; &nbsp; &nbsp;PetscInt &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;l, m, n;//L, M, N<br>&nbsp; } UserContext;<br><br>public:<br>&nbsp; CPETScPoissonSolver(int argc, char** argv);<br>&nbsp; ~CPETScPoissonSolver(void);<br><br>&nbsp; //........<br><br>private:<br><br>&nbsp; //Yale Sparse Matrix format matrix<br>
&nbsp; &nbsp;PetscScalar* &nbsp;A;<br>&nbsp; &nbsp;PetscInt &nbsp; &nbsp; * &nbsp; I;<br>&nbsp; &nbsp;PetscInt &nbsp; &nbsp; * &nbsp; J;<br><br>&nbsp; &nbsp;// Number of Nonzero Element<br>&nbsp; &nbsp;PetscInt &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;nnz;<br>&nbsp; &nbsp;//grid step<br>&nbsp; &nbsp;PetscInt &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;L, M, N;<br>&nbsp; &nbsp;UserContext &nbsp; &nbsp;userctx;<br>
private:<br>&nbsp; bool &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; FirstTime;<br>};<br><br>//in cpp<br>static char helpPetscPoisson[] = &quot;PETSc class Solves a variable Poisson problem with Null Space Method.\n\n&quot;;<br><br>CPETScPoissonSolver::CPETScPoissonSolver(int argc, char** argv)<br>
{<br>&nbsp; PetscInitialize(&amp;argc, &amp;argv, (char*)0, helpPetscPoisson);<br>&nbsp; FirstTime=true;<br>}<br>CPETScPoissonSolver::~CPETScPoissonSolver(void)<br>{<br>&nbsp; PetscFinalize();<br>}<br>......<br>void CPETScPoissonSolver::SetAIJ(PetscScalar *a, PetscInt *i, PetscInt *j, PetscInt Nnz)<br>
{<br>&nbsp; A= a; I=i; J=j; nnz = Nnz;<br>}<br><br>PetscErrorCode CPETScPoissonSolver::UserInitializeLinearSolver()<br>{<br>&nbsp; PetscErrorCode ierr = 0;<br>&nbsp; PetscInt Num = L*M*N;<br><br>&nbsp; //Since we use successive solvers, so in the second time step we must deallocate the original matrix then setup a new one<br>
&nbsp; if(FirstTime==true)<br>&nbsp; {<br>&nbsp; &nbsp; &nbsp; ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, &nbsp;Num, &nbsp;Num, &nbsp;I, &nbsp;J, &nbsp;A, &amp;userctx.A); CHKERRQ(ierr);<br>&nbsp; }<br>&nbsp; else<br>&nbsp; {<br>&nbsp; &nbsp; &nbsp; &nbsp; FirstTime = false;<br>&nbsp; &nbsp; &nbsp; &nbsp; ierr = MatDestroy(userctx.A);CHKERRQ(ierr);<br>
&nbsp; &nbsp; &nbsp; &nbsp; ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, &nbsp;Num, &nbsp;Num, &nbsp;I, &nbsp;J, &nbsp;A, &amp;userctx.A); CHKERRQ(ierr);<br>&nbsp; }<br><br>&nbsp; if(FirstTime==true)<br>&nbsp; {<br><br>&nbsp; &nbsp; &nbsp; ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,Num,PETSC_NULL,&amp;userctx.b);CHKERRQ(ierr);<br>
&nbsp; &nbsp; &nbsp; ierr = VecDuplicate(userctx.b,&amp;userctx.x);CHKERRQ(ierr);<br>&nbsp; &nbsp; &nbsp; ierr = MatNullSpaceCreate(PETSC_COMM_SELF, PETSC_TRUE, 0, PETSC_NULL, &amp;userctx.nullspace); CHKERRQ(ierr);<br>&nbsp; &nbsp; &nbsp; ierr = KSPCreate(PETSC_COMM_SELF,&amp;userctx.ksp);CHKERRQ(ierr);<br>
&nbsp; &nbsp; &nbsp; /*Set Null Space for KSP*/<br>&nbsp; &nbsp; &nbsp; ierr = KSPSetNullSpace(userctx.ksp, userctx.nullspace);CHKERRQ(ierr);<br>&nbsp; }<br>&nbsp; return 0;<br>}<br><br><br>PetscErrorCode CPETScPoissonSolver::UserSetBX(PetscScalar *x, PetscScalar *b)<br>
{<br>&nbsp; &nbsp; PetscErrorCode ierr ;<br>&nbsp; &nbsp; //below code we must set it every time step<br>&nbsp; &nbsp; ierr = VecPlaceArray(userctx.x,x);CHKERRQ(ierr);<br>&nbsp; &nbsp; ierr = VecPlaceArray(userctx.b,b);CHKERRQ(ierr);<br>&nbsp; &nbsp; ierr = &nbsp;MatNullSpaceRemove(userctx.nullspace,userctx.b, PETSC_NULL);CHKERRQ(ierr);<br>
&nbsp; &nbsp; return 0;<br>}<br><br>PetscInt CPETScPoissonSolver::UserSolve()<br>{<br>&nbsp; &nbsp; PetscErrorCode ierr;<br>&nbsp; &nbsp; ierr = KSPSetOperators(userctx.ksp,userctx.A,userctx.A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);<br>&nbsp; &nbsp; ierr = KSPSetType(userctx.ksp, KSPCG);<br>
&nbsp; &nbsp; ierr = KSPSetInitialGuessNonzero(userctx.ksp, PETSC_TRUE);<br>&nbsp; &nbsp; ierr = KSPGetPC(userctx.ksp,&amp;userctx.pc);CHKERRQ(ierr);<br>&nbsp; &nbsp; ierr = PCSetType(userctx.pc,PCICC);CHKERRQ(ierr);<br>&nbsp; &nbsp; ierr = PCFactorSetShiftPd(userctx.pc, PETSC_TRUE);<br>
&nbsp; &nbsp; ierr = KSPSetTolerances(userctx.ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,2000);<br>&nbsp; &nbsp; ierr = KSPSetFromOptions(userctx.ksp);CHKERRQ(ierr);<br>&nbsp; &nbsp; /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -<br>
&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; Solve the linear system<br>&nbsp; &nbsp; &nbsp; &nbsp;- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */<br>&nbsp; &nbsp; ierr = KSPSolve(userctx.ksp,userctx.b,userctx.x);CHKERRQ(ierr);<br><br>&nbsp; &nbsp; ierr = VecResetArray(userctx.x);CHKERRQ(ierr);<br>
&nbsp; &nbsp; ierr = VecResetArray(userctx.b);CHKERRQ(ierr);<br><br>&nbsp; &nbsp; return 0;<br>}<br><br>PetscErrorCode CPETScPoissonSolver::ReleaseMem()<br>{<br>&nbsp; &nbsp; PetscErrorCode ierr;<br>&nbsp; &nbsp; ierr = KSPDestroy(userctx.ksp);CHKERRQ(ierr);<br>
&nbsp; &nbsp; ierr = VecDestroy(userctx.x); &nbsp; &nbsp; CHKERRQ(ierr);<br>&nbsp; &nbsp; ierr = VecDestroy(userctx.b); &nbsp; &nbsp; CHKERRQ(ierr);<br>&nbsp; &nbsp; ierr = MatDestroy(userctx.A); &nbsp; &nbsp;CHKERRQ(ierr);<br>&nbsp; &nbsp; ierr = MatNullSpaceDestroy(userctx.nullspace); CHKERRQ(ierr);<br>
&nbsp; &nbsp; return 0;<br>}<br><br>Thanks very much!<br><br><br>-- <br>Pang Shengyong<br>Solidification Simulation Lab,<br>State Key Lab of Mould &amp; Die Technology,<br>Huazhong Univ. of Sci. &amp; Tech. China<br><br><br><br><br>
-- <br>Pang Shengyong<br>Solidification Simulation Lab,<br>State Key Lab of Mould &amp; Die Technology,<br>Huazhong Univ. of Sci. &amp; Tech. China<br></blockquote><br></div></div></blockquote></div><br><br clear="all"><br>
-- <br>Pang Shengyong<br>Solidification Simulation Lab, <br>State Key Lab of Mould &amp; Die Technology,<br>Huazhong Univ. of Sci. &amp; Tech. China<br>