<div dir="ltr">If do not use PETSc or laspack,the result is fine. I have comparied with some literature results. <br><br>Maybe there are some bugs in my coeffient matrix generation code. I will try to fix it. <br><br>Thanks. <br>
<br><div class="gmail_quote">On Tue, Aug 19, 2008 at 3:23 AM, Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
<br>
Again, what happens with -pc_type lu? Does the simulation run for any number of time steps with correct answers?<br><font color="#888888">
<br>
Barry</font><div><div></div><div class="Wj3C7c"><br>
<br>
On Aug 18, 2008, at 12:02 PM, Shengyong wrote:<br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
Hi, Barry<br>
<br>
Thanks for your help!<br>
<br>
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);<br>
<br>
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.<br>
<br>
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.<br>
<br>
Best Regards,<br>
<br>
<br>
<br>
On Mon, Aug 18, 2008 at 11:24 PM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>> wrote:<br>
<br>
On Aug 18, 2008, at 12:35 AM, Shengyong wrote:<br>
<br>
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>
<br>
What does failed mean? Converged for a while then stagnated? Diverged (the residual norm got larger and larger)? Never converged at all?<br>
<br>
icc(1) also failed.<br>
<br>
What does failed mean?<br>
<br>
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>
<br>
Barry<br>
<br>
<br>
<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. 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 <<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>> wrote:<br>
<br>
Try using -pc_type sor -pc_sor_local_symmetric -ksp_type cg<br>
<br>
Also try running the original icc one with the additional option -ksp_monitor_true_residual, see if funky stuff starts to happen.<br>
<br>
You could also try adding -pc_factor_levels 1 to try ICC(1) instead of ICC(0).<br>
<br>
Your code looks fine, I don't see a problem there,<br>
<br>
Barry<br>
<br>
<br>
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'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. 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.<br>
<br>
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:<br>
<br>
//in Header<br>
class CPETScPoissonSolver<br>
{<br>
typedef struct {<br>
Vec x,b; //x, b<br>
Mat A; //A<br>
KSP ksp; //Krylov subspace preconditioner<br>
PC pc;<br>
MatNullSpace nullspace;<br>
PetscInt l, m, n;//L, M, N<br>
} UserContext;<br>
<br>
public:<br>
CPETScPoissonSolver(int argc, char** argv);<br>
~CPETScPoissonSolver(void);<br>
<br>
//........<br>
<br>
private:<br>
<br>
//Yale Sparse Matrix format matrix<br>
PetscScalar* A;<br>
PetscInt * I;<br>
PetscInt * J;<br>
<br>
// Number of Nonzero Element<br>
PetscInt nnz;<br>
//grid step<br>
PetscInt L, M, N;<br>
UserContext userctx;<br>
private:<br>
bool FirstTime;<br>
};<br>
<br>
//in cpp<br>
static char helpPetscPoisson[] = "PETSc class Solves a variable Poisson problem with Null Space Method.\n\n";<br>
<br>
CPETScPoissonSolver::CPETScPoissonSolver(int argc, char** argv)<br>
{<br>
PetscInitialize(&argc, &argv, (char*)0, helpPetscPoisson);<br>
FirstTime=true;<br>
}<br>
CPETScPoissonSolver::~CPETScPoissonSolver(void)<br>
{<br>
PetscFinalize();<br>
}<br>
......<br>
void CPETScPoissonSolver::SetAIJ(PetscScalar *a, PetscInt *i, PetscInt *j, PetscInt Nnz)<br>
{<br>
A= a; I=i; J=j; nnz = Nnz;<br>
}<br>
<br>
PetscErrorCode CPETScPoissonSolver::UserInitializeLinearSolver()<br>
{<br>
PetscErrorCode ierr = 0;<br>
PetscInt Num = L*M*N;<br>
<br>
//Since we use successive solvers, so in the second time step we must deallocate the original matrix then setup a new one<br>
if(FirstTime==true)<br>
{<br>
ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, Num, Num, I, J, A, &userctx.A); CHKERRQ(ierr);<br>
}<br>
else<br>
{<br>
FirstTime = false;<br>
ierr = MatDestroy(userctx.A);CHKERRQ(ierr);<br>
ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, Num, Num, I, J, A, &userctx.A); CHKERRQ(ierr);<br>
}<br>
<br>
if(FirstTime==true)<br>
{<br>
<br>
ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,Num,PETSC_NULL,&userctx.b);CHKERRQ(ierr);<br>
ierr = VecDuplicate(userctx.b,&userctx.x);CHKERRQ(ierr);<br>
ierr = MatNullSpaceCreate(PETSC_COMM_SELF, PETSC_TRUE, 0, PETSC_NULL, &userctx.nullspace); CHKERRQ(ierr);<br>
ierr = KSPCreate(PETSC_COMM_SELF,&userctx.ksp);CHKERRQ(ierr);<br>
/*Set Null Space for KSP*/<br>
ierr = KSPSetNullSpace(userctx.ksp, userctx.nullspace);CHKERRQ(ierr);<br>
}<br>
return 0;<br>
}<br>
<br>
<br>
PetscErrorCode CPETScPoissonSolver::UserSetBX(PetscScalar *x, PetscScalar *b)<br>
{<br>
PetscErrorCode ierr ;<br>
//below code we must set it every time step<br>
ierr = VecPlaceArray(userctx.x,x);CHKERRQ(ierr);<br>
ierr = VecPlaceArray(userctx.b,b);CHKERRQ(ierr);<br>
ierr = MatNullSpaceRemove(userctx.nullspace,userctx.b, PETSC_NULL);CHKERRQ(ierr);<br>
return 0;<br>
}<br>
<br>
PetscInt CPETScPoissonSolver::UserSolve()<br>
{<br>
PetscErrorCode ierr;<br>
ierr = KSPSetOperators(userctx.ksp,userctx.A,userctx.A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);<br>
ierr = KSPSetType(userctx.ksp, KSPCG);<br>
ierr = KSPSetInitialGuessNonzero(userctx.ksp, PETSC_TRUE);<br>
ierr = KSPGetPC(userctx.ksp,&userctx.pc);CHKERRQ(ierr);<br>
ierr = PCSetType(userctx.pc,PCICC);CHKERRQ(ierr);<br>
ierr = PCFactorSetShiftPd(userctx.pc, PETSC_TRUE);<br>
ierr = KSPSetTolerances(userctx.ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,2000);<br>
ierr = KSPSetFromOptions(userctx.ksp);CHKERRQ(ierr);<br>
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -<br>
Solve the linear system<br>
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */<br>
ierr = KSPSolve(userctx.ksp,userctx.b,userctx.x);CHKERRQ(ierr);<br>
<br>
ierr = VecResetArray(userctx.x);CHKERRQ(ierr);<br>
ierr = VecResetArray(userctx.b);CHKERRQ(ierr);<br>
<br>
return 0;<br>
}<br>
<br>
PetscErrorCode CPETScPoissonSolver::ReleaseMem()<br>
{<br>
PetscErrorCode ierr;<br>
ierr = KSPDestroy(userctx.ksp);CHKERRQ(ierr);<br>
ierr = VecDestroy(userctx.x); CHKERRQ(ierr);<br>
ierr = VecDestroy(userctx.b); CHKERRQ(ierr);<br>
ierr = MatDestroy(userctx.A); CHKERRQ(ierr);<br>
ierr = MatNullSpaceDestroy(userctx.nullspace); CHKERRQ(ierr);<br>
return 0;<br>
}<br>
<br>
Thanks very much!<br>
<br>
<br>
-- <br>
Pang Shengyong<br>
Solidification Simulation Lab,<br>
State Key Lab of Mould & Die Technology,<br>
Huazhong Univ. of Sci. & Tech. China<br>
<br>
<br>
<br>
<br>
-- <br>
Pang Shengyong<br>
Solidification Simulation Lab,<br>
State Key Lab of Mould & Die Technology,<br>
Huazhong Univ. of Sci. & Tech. China<br>
<br>
<br>
<br>
<br>
-- <br>
Pang Shengyong<br>
Solidification Simulation Lab,<br>
State Key Lab of Mould & Die Technology,<br>
Huazhong Univ. of Sci. & 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 & Die Technology,<br>Huazhong Univ. of Sci. & Tech. China<br>
</div>