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;}