Actual source code: ex47.c

  1: static char help[] = "Solves -Laplacian u - exp(u) = 0,  0 < x < 1\n\n";
 2:  #include petscda.h
 3:  #include petscsnes.h

  6: int main(int argc,char **argv) {
  7:   SNES snes; Vec x,f; Mat J; DA da;
  8:   PetscInitialize(&argc,&argv,(char *)0,help);

 10:   DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,8,1,1,PETSC_NULL,&da);
 11:   DACreateGlobalVector(da,&x); VecDuplicate(x,&f);
 12:   DAGetMatrix(da,MATAIJ,&J);

 14:   SNESCreate(PETSC_COMM_WORLD,&snes);
 15:   SNESSetFunction(snes,f,ComputeFunction,da);
 16:   SNESSetJacobian(snes,J,J,ComputeJacobian,da);
 17:   SNESSetFromOptions(snes);
 18:   SNESSolve(snes,PETSC_NULL,x);

 20:   MatDestroy(J); VecDestroy(x); VecDestroy(f); SNESDestroy(snes); DADestroy(da);
 21:   PetscFinalize();
 22:   return 0;}
 23: PetscErrorCode ComputeFunction(SNES snes,Vec x,Vec f,void *ctx) {
 24:   PetscInt i,Mx,xs,xm; PetscScalar *xx,*ff,hx; DA da = (DA) ctx; Vec xlocal;
 25:   DAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
 26:   hx     = 1.0/(PetscReal)(Mx-1);
 27:   DAGetLocalVector(da,&xlocal);DAGlobalToLocalBegin(da,x,INSERT_VALUES,xlocal);  DAGlobalToLocalEnd(da,x,INSERT_VALUES,xlocal);
 28:   DAVecGetArray(da,xlocal,&xx); DAVecGetArray(da,f,&ff);
 29:   DAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);

 31:   for (i=xs; i<xs+xm; i++) {
 32:     if (i == 0 || i == Mx-1) ff[i] = xx[i]/hx;
 33:     else  ff[i] =  (2.0*xx[i] - xx[i-1] - xx[i+1])/hx - hx*PetscExpScalar(xx[i]);
 34:   }
 35:   DAVecRestoreArray(da,xlocal,&xx); DARestoreLocalVector(da,&xlocal);DAVecRestoreArray(da,f,&ff);
 36:   return 0;}
 37: PetscErrorCode ComputeJacobian(SNES snes,Vec x,Mat *J,Mat *B,MatStructure *flag,void *ctx){
 38:   DA da = (DA) ctx; PetscInt i,Mx,xm,xs; PetscScalar hx,*xx; Vec xlocal;
 39:   DAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
 40:   hx = 1.0/(PetscReal)(Mx-1);
 41:   DAGetLocalVector(da,&xlocal);DAGlobalToLocalBegin(da,x,INSERT_VALUES,xlocal);  DAGlobalToLocalEnd(da,x,INSERT_VALUES,xlocal);
 42:   DAVecGetArray(da,xlocal,&xx);
 43:   DAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);

 45:   for (i=xs; i<xs+xm; i++) {
 46:     if (i == 0 || i == Mx-1) { MatSetValue(*J,i,i,1.0/hx,INSERT_VALUES);}
 47:     else {
 48:       MatSetValue(*J,i,i-1,-1.0/hx,INSERT_VALUES);
 49:       MatSetValue(*J,i,i,2.0/hx - hx*PetscExpScalar(xx[i]),INSERT_VALUES);
 50:       MatSetValue(*J,i,i+1,-1.0/hx,INSERT_VALUES);
 51:     }
 52:   }
 53:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY); MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);  *flag = SAME_NONZERO_PATTERN;
 54:   DAVecRestoreArray(da,xlocal,&xx);DARestoreLocalVector(da,&xlocal);
 55:   return 0;}