[petsc-users] Problem in SNES
Ali Reza Khaz'ali
arkhazali at cc.iut.ac.ir
Thu Jul 12 10:14:06 CDT 2018
Thanks very much for your answer,
I am experimenting with SNES for now, however, my purpose is to solve a
multi-physics large system with unstructured grids (Enhanced Oil
Recovery simulation in fractured petroleum reservoirs to be specific,
phase equilibrium and three phase flow equations must be solved
simultaneously), where to the extend of my knowledge, no existing PETSc
vector structure can be used. Therefore, I started with solving simple
ODEs with SNES without using those structures.
An interesting fact about my code samples is, with increasing the number
of MPI workers, the number of SNES iterations (and hence, the precision
of the solution) decrease. I'm sure there is a flaw in my code (I do not
have much MPI experience), but I can't find it.
Best wishes,
Ali
On 7/12/2018 3:58 PM, Matthew Knepley wrote:
> On Thu, Jul 12, 2018 at 6:30 AM Ali Reza Khaz'ali
> <arkhazali at cc.iut.ac.ir <mailto:arkhazali at cc.iut.ac.ir>> wrote:
>
> Hi,
>
> I have to apologize if I'm asking a noob question. I am trying to
> solve
> y"-(1-x/5)yy'=x with boundary condition of y(1)=2 and y(3)=-1 with
> finite difference method. The discretized equations is to be
> solved by
> SNES. I wrote the follwoing code, however, it only works if number of
> MPI processes is set to one. for more workers, it does not seem to
> give
> valid results. I'd be grateful if anybody can guide me to debug the
> code, please.
>
>
> Why not just use DMDA? There are lots of examples (SNES ex5 for
> instance) and it
> handles the parallelism for you. This problem could be done in a few
> hours I think.
>
> Thanks,
>
> Matt
>
> Regards,
>
> Ali
>
>
>
> static char help[] = "Newton's method for a non-linear ODE
> system.\n\n";
> double h;
> int nParts;
>
> #include <petscsnes.h>
>
> /*
> User-defined routines
> */
> extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void*);
> extern PetscErrorCode FormFunction(SNES, Vec, Vec, void*);
>
> int main(int argc, char **argv)
> {
> SNES snes; /* nonlinear solver context */
> KSP ksp; /* linear solver context */
> PC pc; /* preconditioner context */
> Vec x, r; /* solution, residual vectors */
> Mat J; /* Jacobian matrix */
> PetscErrorCode ierr;
> PetscMPIInt size;
> PetscScalar pfive = 0, *xx;
> int rank;
> int i, j;
> Vec f;
>
>
> nParts = 30;
> h = (3.0 - 1.0) / nParts;
>
> ierr = PetscInitialize(&argc, &argv, (char*)0, help); if (ierr)
> return ierr;
> ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRQ(ierr);
> ierr = SNESCreate(PETSC_COMM_WORLD, &snes); CHKERRQ(ierr);
> ierr = VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, nParts-1,
> &x);
> CHKERRQ(ierr);
> ierr = VecDuplicate(x, &r); CHKERRQ(ierr);
> ierr = MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE,
> PETSC_DECIDE,
> nParts - 1, nParts - 1, 3, NULL, 3, NULL, &J); CHKERRQ(ierr);
>
>
>
>
> ierr = SNESSetFunction(snes, r, FormFunction, NULL);
> CHKERRQ(ierr);
> ierr = SNESSetJacobian(snes, J, J, FormJacobian, NULL);
> CHKERRQ(ierr);
>
>
>
> ierr = SNESGetKSP(snes, &ksp); CHKERRQ(ierr);
> ierr = KSPGetPC(ksp, &pc); CHKERRQ(ierr);
> ierr = KSPSetTolerances(ksp, 1.e-8, 1e-16, PETSC_DEFAULT, 200);
> CHKERRQ(ierr);
> ierr = SNESSetTolerances(snes, 1.e-8, 1e-16, PETSC_DEFAULT, 200,
> 200); CHKERRQ(ierr);
> ierr = SNESSetType(snes, SNESNEWTONLS); CHKERRQ(ierr);
> ierr = KSPSetType(ksp, KSPCG); CHKERRQ(ierr);
>
> ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);
> ierr = VecSet(x, pfive); CHKERRQ(ierr);
> ierr = SNESSolve(snes, NULL, x); CHKERRQ(ierr);
> ierr = SNESGetFunction(snes, &f, 0, 0); CHKERRQ(ierr);
>
>
>
> ierr = VecDestroy(&x); CHKERRQ(ierr); ierr = VecDestroy(&r);
> CHKERRQ(ierr);
> ierr = MatDestroy(&J); CHKERRQ(ierr); ierr = SNESDestroy(&snes);
> CHKERRQ(ierr);
> ierr = PetscFinalize();
> return ierr;
> }
>
> PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *ctx)
> {
> PetscErrorCode ierr;
> const PetscScalar *xx;
> PetscScalar *ff;
> PetscInt vstart, vend;
> int i, j, rank, size, xlength;
> double xp, xxstart, xxend, x1, x2;
> MPI_Status MyStat;
>
> xxstart = 0;
> xxend = 0;
>
>
> MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
> MPI_Comm_size(PETSC_COMM_WORLD, &size);
>
> ierr = VecGetOwnershipRange(x, &vstart, &vend); CHKERRQ(ierr);
> xlength = vend - vstart;
> ierr = VecGetArrayRead(x, &xx); CHKERRQ(ierr);
> ierr = VecGetArray(f, &ff); CHKERRQ(ierr);
> ierr = VecView(x, PETSC_VIEWER_STDOUT_WORLD);
> if (size > 1) {
> if (rank == 0) {
> ierr = MPI_Send(&xx[xlength - 1], 1, MPI_DOUBLE, 1, 0,
> PETSC_COMM_WORLD);
> ierr = MPI_Recv(&xxend, 1, MPI_DOUBLE, 1, 1,
> PETSC_COMM_WORLD, &MyStat);
> }
> else if (rank == (size - 1)) {
> ierr = MPI_Recv(&xxstart, 1, MPI_DOUBLE, rank - 1,
> rank -
> 1, PETSC_COMM_WORLD, &MyStat);
> ierr = MPI_Send(&xx[0], 1, MPI_DOUBLE, rank - 1, rank,
> PETSC_COMM_WORLD);
> }
> else {
> ierr = MPI_Send(&xx[xlength - 1], 1, MPI_DOUBLE, rank
> + 1,
> rank, PETSC_COMM_WORLD);
> ierr = MPI_Recv(&xxend, 1, MPI_DOUBLE, rank + 1, rank
> + 1,
> PETSC_COMM_WORLD, &MyStat);
> ierr = MPI_Send(&xx[0], 1, MPI_DOUBLE, rank - 1, rank,
> PETSC_COMM_WORLD);
> ierr = MPI_Recv(&xxstart, 1, MPI_DOUBLE, rank - 1,
> rank -
> 1, PETSC_COMM_WORLD, &MyStat);
> }
> }
>
>
>
>
> for (i = vstart; i < vend; i++) {
> xp = 1 + i * h;
> j = i - vstart;
>
>
> if (i == vstart) x1 = xxstart;
> else x1 = xx[j - 1];
> if (i == (vend - 1)) x2 = xxend;
> else x2 = xx[j + 1];
>
>
> if (i==0) ff[j] = (x2 - 2 * xx[j] + 2) / (h*h) - (1 - xp /
> 5)*xx[j] * (x2 - 2) / (2 * h) - xp;
> else if (i==(nParts-2)) ff[j]=(-1 - 2 * xx[j] + x1) /
> (h*h) -
> (1 - xp / 5)*xx[j] * (-1 - x1) / (2 * h) - xp;
> else ff[j] = (x2 - 2 * xx[j] + x1) / (h*h) - (1 - xp /
> 5)*xx[j]
> * (x2 - x1) / (2 * h) - xp;
> }
>
> /* Restore vectors */
> ierr = VecRestoreArrayRead(x, &xx); CHKERRQ(ierr);
> ierr = VecRestoreArray(f, &ff); CHKERRQ(ierr);
> return 0;
> }
>
> PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void
> *dummy)
> {
> const PetscScalar *xx;
> PetscErrorCode ierr;
>
> PetscInt vstart, vend;
> double rrow[3];
> int col[3];
> int i, j, rank, size, xlength;
> double xp, xxstart, xxend, x1, x2;
> MPI_Status MyStat;
>
> ierr = MatZeroEntries(jac); CHKERRQ(ierr);
>
> ierr = VecGetOwnershipRange(x, &vstart, &vend); CHKERRQ(ierr);
> xlength = vend - vstart;
> MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
> MPI_Comm_size(PETSC_COMM_WORLD, &size);
>
> xxstart = 0;
> xxend = 0;
>
> /*
> Get pointer to vector data.0+1
> */
> ierr = VecGetArrayRead(x, &xx); CHKERRQ(ierr);
> if (size > 1) {
> if (rank == 0) {
> ierr = MPI_Send(&xx[xlength - 1], 1, MPI_DOUBLE, 1, 0,
> PETSC_COMM_WORLD);
> ierr = MPI_Recv(&xxend, 1, MPI_DOUBLE, 1, 1,
> PETSC_COMM_WORLD, &MyStat);
> xxstart = 1e12;
> }
> else if (rank == (size - 1)) {
> ierr = MPI_Recv(&xxstart, 1, MPI_DOUBLE, rank - 1,
> rank -
> 1, PETSC_COMM_WORLD, &MyStat);
> ierr = MPI_Send(&xx[0], 1, MPI_DOUBLE, rank - 1, rank,
> PETSC_COMM_WORLD);
> xxend = 1e12;
> }
> else {
> ierr = MPI_Send(&xx[xlength - 1], 1, MPI_DOUBLE, rank
> + 1,
> rank, PETSC_COMM_WORLD);
> ierr = MPI_Recv(&xxend, 1, MPI_DOUBLE, rank + 1, rank
> + 1,
> PETSC_COMM_WORLD, &MyStat);
> ierr = MPI_Send(&xx[0], 1, MPI_DOUBLE, rank - 1, rank,
> PETSC_COMM_WORLD);
> ierr = MPI_Recv(&xxstart, 1, MPI_DOUBLE, rank - 1,
> rank -
> 1, PETSC_COMM_WORLD, &MyStat);
> }
> }
>
>
> for (i = vstart; i < vend; i++) {
> xp = 1 + i * h;
> j = i - vstart;
>
> if (i == vstart) x1 = xxstart;
> else x1 = xx[j - 1];
> if (i == (vend - 1)) x2 = xxend;
> else x2 = xx[j + 1];
>
> if (i == 0) {
> rrow[0] = -2 / (h*h) - (1 - xp / 5)*(x2 - 2) / (2 * h);
> rrow[1] = 1 / (h*h) - (1 + xp / 5)*xx[j] / (2 * h);
> col[0] = i;
> col[1] = i + 1;
> ierr = MatSetValues(jac, 1, &i, 2, col, rrow,
> INSERT_VALUES); CHKERRQ(ierr);
> }
> else if (i == (nParts-2)) {
> rrow[0] = 1 / (h*h) + (1 + xp / 5)*xx[j] / (2 * h);
> rrow[1] = -2 / (h*h) - (1 - xp / 5)*(-1 - x1) / (2 * h);
> col[0] = i - 1;
> col[1] = i;
> ierr = MatSetValues(jac, 1, &i, 2, col, rrow,
> INSERT_VALUES); CHKERRQ(ierr);
> }
> else {
> rrow[0] = 1 / (h*h) + (1 + xp / 5)*xx[j] / (2 * h);
> rrow[1] = -2 / (h*h) - (1 - xp / 5)*(x2 - x1) / (2 * h);
> rrow[2] = 1 / (h*h) - (1 + xp / 5)*xx[j] / (2 * h);
> col[0] = i-1;
> col[1] = i;
> col[2] = i + 1;
> ierr = MatSetValues(jac, 1, &i, 3, col, rrow,
> INSERT_VALUES); CHKERRQ(ierr);
> }
> }
>
>
>
> /*
> Restore vector
> */
> ierr = VecRestoreArrayRead(x, &xx); CHKERRQ(ierr);
>
> /*
> Assemble matrix
> */
> B = jac;
> ierr = MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
> ierr = MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
> ierr = MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
> ierr = MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
>
> return 0;
> }
>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which
> their experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/ <http://www.caam.rice.edu/%7Emk51/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180712/9da7481b/attachment-0001.html>
More information about the petsc-users
mailing list