<div dir="ltr">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 clear="all"><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>
</div>