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: [email protected] 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 <[email protected]> wrote: > On Thu, Jul 12, 2018 at 11:14 AM Ali Reza Khaz'ali <[email protected]> > 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 <[email protected]> >> 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/> >
