static char help[] = "Bratu nonlinear PDE in 3d.\n\ We solve the Bratu (SFI - solid fuel ignition) problem in a 3D rectangular\n\ domain, using distributed arrays (DMDAs) to partition the parallel grid.\n\ The command line options include:\n\ -par , where indicates the problem's nonlinearity\n\ problem SFI: = Bratu parameter (0 <= par <= 6.81)\n\n"; /*T Concepts: SNES^parallel Bratu example Concepts: DMDA^using distributed arrays; Processors: n T*/ /* ------------------------------------------------------------------------ Solid Fuel Ignition (SFI) problem. This problem is modeled by the partial differential equation -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1, with boundary conditions u = 0 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1 A finite difference approximation with the usual 7-point stencil is used to discretize the boundary value problem to obtain a nonlinear system of equations. ------------------------------------------------------------------------- */ /* Include "petscdmda.h" so that we can use distributed arrays (DMDAs). Include "petscsnes.h" so that we can use SNES solvers. Note that this file automatically includes: petscsys.h - base PETSc routines petscvec.h - vectors petscmat.h - matrices petscis.h - index sets petscksp.h - Krylov subspace methods petscviewer.h - viewers petscpc.h - preconditioners petscksp.h - linear solvers */ #include #include #include /* User-defined application context - contains data needed by the application-provided call-back routines, FormJacobian() and FormFunction(). */ typedef struct { PetscReal param; /* test problem parameter */ DM da; /* distributed array data structure */ } AppCtx; /* User-defined routines */ extern PetscErrorCode FormFunctionLocal(SNES,Vec,Vec,void*); extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*); extern PetscErrorCode FormInitialGuess(AppCtx*,Vec); extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*); int main(int argc,char **argv) { SNES snes; /* nonlinear solver */ Vec x,r; /* solution, residual vectors */ Mat J = NULL; /* Jacobian matrix */ AppCtx user; /* user-defined work context */ PetscInt its; /* iterations for convergence */ MatFDColoring matfdcoloring = NULL; PetscBool matrix_free = PETSC_FALSE,coloring = PETSC_FALSE, coloring_ds = PETSC_FALSE,local_coloring = PETSC_FALSE; PetscErrorCode ierr; PetscReal bratu_lambda_max = 6.81,bratu_lambda_min = 0.,fnorm; /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Initialize program - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Initialize problem parameters - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ user.param = 6.0; ierr = PetscOptionsGetReal(NULL,NULL,"-par",&user.param,NULL);CHKERRQ(ierr); if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) SETERRQ(PETSC_COMM_SELF,1,"Lambda is out of range"); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Create nonlinear solver context - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Create distributed array (DMDA) to manage parallel grid and vectors - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,4,4,4,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,NULL,&user.da);CHKERRQ(ierr); ierr = DMSetFromOptions(user.da);CHKERRQ(ierr); ierr = DMSetUp(user.da);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Extract global vectors from DMDA; then duplicate for remaining vectors that are the same types - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = DMCreateGlobalVector(user.da,&x);CHKERRQ(ierr); ierr = VecDuplicate(x,&r);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Set function evaluation routine and vector - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = SNESSetFunction(snes,r,FormFunction,(void*)&user);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Create matrix data structure; set Jacobian evaluation routine Set Jacobian matrix data structure and default Jacobian evaluation routine. User can override with: -snes_mf : matrix-free Newton-Krylov method with no preconditioning (unless user explicitly sets preconditioner) -snes_mf_operator : form preconditioning matrix as set by the user, but use matrix-free approx for Jacobian-vector products within Newton-Krylov method -fdcoloring : using finite differences with coloring to compute the Jacobian Note one can use -matfd_coloring wp or ds the only reason for the -fdcoloring_ds option below is to test the call to MatFDColoringSetType(). - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = PetscOptionsGetBool(NULL,NULL,"-snes_mf",&matrix_free,NULL);CHKERRQ(ierr); ierr = PetscOptionsGetBool(NULL,NULL,"-fdcoloring",&coloring,NULL);CHKERRQ(ierr); ierr = PetscOptionsGetBool(NULL,NULL,"-fdcoloring_ds",&coloring_ds,NULL);CHKERRQ(ierr); ierr = PetscOptionsGetBool(NULL,NULL,"-fdcoloring_local",&local_coloring,NULL);CHKERRQ(ierr); if (!matrix_free) { ierr = DMSetMatType(user.da,MATAIJ);CHKERRQ(ierr); ierr = DMCreateMatrix(user.da,&J);CHKERRQ(ierr); if (coloring) { ISColoring iscoloring; if (!local_coloring) { ierr = DMCreateColoring(user.da,IS_COLORING_GLOBAL,&iscoloring);CHKERRQ(ierr); ierr = MatFDColoringCreate(J,iscoloring,&matfdcoloring);CHKERRQ(ierr); ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunction,&user);CHKERRQ(ierr); } else { ierr = DMCreateColoring(user.da,IS_COLORING_LOCAL,&iscoloring);CHKERRQ(ierr); ierr = MatFDColoringCreate(J,iscoloring,&matfdcoloring);CHKERRQ(ierr); ierr = MatFDColoringUseDM(J,matfdcoloring);CHKERRQ(ierr); ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunctionLocal,&user);CHKERRQ(ierr); } if (coloring_ds) { ierr = MatFDColoringSetType(matfdcoloring,MATMFFD_DS);CHKERRQ(ierr); } ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); ierr = MatFDColoringSetUp(J,iscoloring,matfdcoloring);CHKERRQ(ierr); ierr = SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,matfdcoloring);CHKERRQ(ierr); ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); } else { ierr = SNESSetJacobian(snes,J,J,FormJacobian,&user);CHKERRQ(ierr); } } /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Customize nonlinear solver; set runtime options - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = SNESSetDM(snes,user.da);CHKERRQ(ierr); ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Evaluate initial guess Note: The user should initialize the vector, x, with the initial guess for the nonlinear solver prior to calling SNESSolve(). In particular, to employ an initial guess of zero, the user should explicitly set this vector to zero by calling VecSet(). */ ierr = FormInitialGuess(&user,x);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Solve nonlinear system - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = SNESSolve(snes,NULL,x);CHKERRQ(ierr); ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Explicitly check norm of the residual of the solution - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = FormFunction(snes,x,r,(void*)&user);CHKERRQ(ierr); ierr = VecNorm(r,NORM_2,&fnorm);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D fnorm %g\n",its,(double)fnorm);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Free work space. All PETSc objects should be destroyed when they are no longer needed. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = MatDestroy(&J);CHKERRQ(ierr); ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&r);CHKERRQ(ierr); ierr = SNESDestroy(&snes);CHKERRQ(ierr); ierr = DMDestroy(&user.da);CHKERRQ(ierr); ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr); ierr = PetscFinalize(); return ierr; } /* ------------------------------------------------------------------- */ /* FormInitialGuess - Forms initial approximation. Input Parameters: user - user-defined application context X - vector Output Parameter: X - vector */ PetscErrorCode FormInitialGuess(AppCtx *user,Vec X) { PetscInt i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm; PetscErrorCode ierr; PetscReal lambda,temp1,hx,hy,hz,tempk,tempj; PetscScalar ***x; PetscFunctionBeginUser; ierr = DMDAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,&Mz,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr); lambda = user->param; hx = 1.0/(PetscReal)(Mx-1); hy = 1.0/(PetscReal)(My-1); hz = 1.0/(PetscReal)(Mz-1); temp1 = lambda/(lambda + 1.0); /* Get a pointer to vector data. - For default PETSc vectors, VecGetArray() returns a pointer to the data array. Otherwise, the routine is implementation dependent. - You MUST call VecRestoreArray() when you no longer need access to the array. */ ierr = DMDAVecGetArray(user->da,X,&x);CHKERRQ(ierr); /* Get local grid boundaries (for 3-dimensional DMDA): xs, ys, zs - starting grid indices (no ghost points) xm, ym, zm - widths of local grid (no ghost points) */ ierr = DMDAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr); /* Compute initial guess over the locally owned part of the grid */ for (k=zs; kda,X,&x);CHKERRQ(ierr); PetscFunctionReturn(0); } /* ------------------------------------------------------------------- */ /* FormFunctionLocal - Evaluates nonlinear function, F(x) on a ghosted patch Input Parameters: . snes - the SNES context . localX - input vector, this contains the ghosted region needed . ptr - optional user-defined context, as set by SNESSetFunction() Output Parameter: . F - function vector, this does not contain a ghosted region */ PetscErrorCode FormFunctionLocal(SNES snes,Vec localX,Vec F,void *ptr) { AppCtx *user = (AppCtx*)ptr; PetscErrorCode ierr; PetscInt i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm; PetscReal two = 2.0,lambda,hx,hy,hz,hxhzdhy,hyhzdhx,hxhydhz,sc; PetscScalar u_north,u_south,u_east,u_west,u_up,u_down,u,u_xx,u_yy,u_zz,***x,***f; DM da; PetscFunctionBeginUser; ierr = SNESGetDM(snes,&da);CHKERRQ(ierr); ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,&Mz,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr); lambda = user->param; hx = 1.0/(PetscReal)(Mx-1); hy = 1.0/(PetscReal)(My-1); hz = 1.0/(PetscReal)(Mz-1); sc = hx*hy*hz*lambda; hxhzdhy = hx*hz/hy; hyhzdhx = hy*hz/hx; hxhydhz = hx*hy/hz; /* Get pointers to vector data */ ierr = DMDAVecGetArrayRead(da,localX,&x);CHKERRQ(ierr); ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); /* Get local grid boundaries */ ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr); /* Compute function over the locally owned part of the grid */ for (k=zs; kparam; hx = 1.0/(PetscReal)(Mx-1); hy = 1.0/(PetscReal)(My-1); hz = 1.0/(PetscReal)(Mz-1); sc = hx*hy*hz*lambda; hxhzdhy = hx*hz/hy; hyhzdhx = hy*hz/hx; hxhydhz = hx*hy/hz; /* Scatter ghost points to local vector, using the 2-step process DMGlobalToLocalBegin(), DMGlobalToLocalEnd(). By placing code between these two statements, computations can be done while messages are in transition. */ ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); /* Get pointer to vector data */ ierr = DMDAVecGetArrayRead(da,localX,&x);CHKERRQ(ierr); /* Get local grid boundaries */ ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr); /* Compute entries for the locally owned part of the Jacobian. - Currently, all PETSc parallel matrix formats are partitioned by contiguous chunks of rows across the processors. - Each processor needs to insert only elements that it owns locally (but any non-local elements will be sent to the appropriate processor during matrix assembly). - Here, we set all entries for a particular row at once. - We can set matrix entries either using either MatSetValuesLocal() or MatSetValues(), as discussed above. */ for (k=zs; k