[petsc-users] Problem in SNES
Adam Denchfield
adenchfi at hawk.iit.edu
Thu Jul 12 11:12:03 CDT 2018
You may want to look at DMNetwork since you are doing a multi-physics
system, http://www.mcs.anl.gov/papers/P7065-0617.pdf.
Regards,
*Adam Denchfield*
*Peer Career Coach - Career Services*
Illinois Institute of Technology
*Bachelors of Science in Applied Physics (2018)*
Email: adenchfi at hawk.iit.edu
My LinkedIn <http://www.linkedin.com/in/adamrobertdenchfield> My
ResearchGate Profile <https://www.researchgate.net/profile/Adam_Denchfield>
On Thu, Jul 12, 2018 at 10:30 AM, Matthew Knepley <knepley at gmail.com> wrote:
> On Thu, Jul 12, 2018 at 11:14 AM Ali Reza Khaz'ali <arkhazali at cc.iut.ac.ir>
> wrote:
>
>> 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.
>>
> The intention is that you use DMPlex for this.
>
>> 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.
>>
> We really do not have the resources to debug raw MPI code. That is why we
> have the examples.
>
> Thanks,
>
> Matt
>
>> 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>
>> 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/>
>>
>>
>>
>
> --
> 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/~mk51/>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180712/e4b2848e/attachment-0001.html>
More information about the petsc-users
mailing list