Actual source code: ex25.c

      2: static char help[] ="Minimum surface problem\n\
      3: Uses 2-dimensional distributed arrays.\n\
      4: \n\
      5:   Solves the linear systems via multilevel methods \n\
      6: \n\n";

      8: /*T
      9:    Concepts: SNES^solving a system of nonlinear equations
     10:    Concepts: DA^using distributed arrays
     11:    Concepts: multigrid;
     12:    Processors: n
     13: T*/

     15: /*
     16:
     17:     This example models the partial differential equation
     18:
     19:          - Div((1 + ||GRAD T||^2)^(1/2) (GRAD T)) = 0.
     20:
     21:
     22:     in the unit square, which is uniformly discretized in each of x and
     23:     y in this simple encoding.  The degrees of freedom are vertex centered
     24:
     25:     A finite difference approximation with the usual 5-point stencil
     26:     is used to discretize the boundary value problem to obtain a
     27:     nonlinear system of equations.
     28:
     29: */

     31:  #include petscsnes.h
     32:  #include petscda.h
     33:  #include petscmg.h
     34:  #include petscdmmg.h

     41: int main(int argc,char **argv)
     42: {
     43:   DMMG           *dmmg;
     44:   SNES           snes;
     46:   PetscInt       its,lits;
     47:   PetscReal      litspit;
     48:   DA             da;

     50:   PetscInitialize(&argc,&argv,PETSC_NULL,help);

     53:   /*
     54:       Create the multilevel DA data structure
     55:   */
     56:   DMMGCreate(PETSC_COMM_WORLD,3,0,&dmmg);

     58:   /*
     59:       Set the DA (grid structure) for the grids.
     60:   */
     61:   DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,-5,-5,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
     62:   DMMGSetDM(dmmg,(DM)da);
     63:   DADestroy(da);

     65:   /*
     66:        Process adiC(36): FormFunctionLocal FormFunctionLocali

     68:      Create the nonlinear solver, and tell the DMMG structure to use it
     69:   */
     70:   /*  DMMGSetSNES(dmmg,FormFunction,0); */
     71:   DMMGSetSNESLocal(dmmg,FormFunctionLocal,0,ad_FormFunctionLocal,0);

     73:   /*
     74:       PreLoadBegin() means that the following section of code is run twice. The first time
     75:      through the flag PreLoading is on this the nonlinear solver is only run for a single step.
     76:      The second time through (the actually timed code) the maximum iterations is set to 10
     77:      Preload of the executable is done to eliminate from the timing the time spent bring the
     78:      executable into memory from disk (paging in).
     79:   */
     80:   PreLoadBegin(PETSC_TRUE,"Solve");
     81:     DMMGSolve(dmmg);
     82:   PreLoadEnd();
     83:   snes = DMMGGetSNES(dmmg);
     84:   SNESGetIterationNumber(snes,&its);
     85:   SNESGetNumberLinearIterations(snes,&lits);
     86:   litspit = ((PetscReal)lits)/((PetscReal)its);
     87:   PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D\n",its);
     88:   PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);
     89:   PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / Newton = %e\n",litspit);

     91:   DMMGDestroy(dmmg);
     92:   PetscFinalize();

     94:   return 0;
     95: }
     96: /* --------------------  Evaluate Function F(x) --------------------- */
     99: PetscErrorCode FormFunction(SNES snes,Vec T,Vec F,void* ptr)
    100: {
    101:   DMMG           dmmg = (DMMG)ptr;
    103:   PetscInt       i,j,mx,my,xs,ys,xm,ym;
    104:   PetscScalar    hx,hy;
    105:   PetscScalar    **t,**f,gradup,graddown,gradleft,gradright,gradx,grady;
    106:   PetscScalar    coeffup,coeffdown,coeffleft,coeffright;
    107:   Vec            localT;

    110:   DAGetLocalVector((DA)dmmg->dm,&localT);
    111:   DAGetInfo((DA)dmmg->dm,PETSC_NULL,&mx,&my,0,0,0,0,0,0,0,0);
    112:   hx    = 1.0/(PetscReal)(mx-1);  hy    = 1.0/(PetscReal)(my-1);
    113:
    114:   /* Get ghost points */
    115:   DAGlobalToLocalBegin((DA)dmmg->dm,T,INSERT_VALUES,localT);
    116:   DAGlobalToLocalEnd((DA)dmmg->dm,T,INSERT_VALUES,localT);
    117:   DAGetCorners((DA)dmmg->dm,&xs,&ys,0,&xm,&ym,0);
    118:   DAVecGetArray((DA)dmmg->dm,localT,&t);
    119:   DAVecGetArray((DA)dmmg->dm,F,&f);

    121:   /* Evaluate function */
    122:   for (j=ys; j<ys+ym; j++) {
    123:     for (i=xs; i<xs+xm; i++) {

    125:       if (i == 0 || i == mx-1 || j == 0 || j == my-1) {

    127:         f[j][i] = t[j][i] - (1.0 - (2.0*hx*(PetscReal)i - 1.0)*(2.0*hx*(PetscReal)i - 1.0));
    128:
    129:       } else {

    131:         gradup     = (t[j+1][i] - t[j][i])/hy;
    132:         graddown   = (t[j][i] - t[j-1][i])/hy;
    133:         gradright  = (t[j][i+1] - t[j][i])/hx;
    134:         gradleft   = (t[j][i] - t[j][i-1])/hx;

    136:         gradx      = .5*(t[j][i+1] - t[j][i-1])/hx;
    137:         grady      = .5*(t[j+1][i] - t[j-1][i])/hy;

    139:         coeffup    = 1.0/PetscSqrtScalar(1.0 + gradup*gradup + gradx*gradx);
    140:         coeffdown  = 1.0/PetscSqrtScalar(1.0 + graddown*graddown + gradx*gradx);

    142:         coeffleft  = 1.0/PetscSqrtScalar(1.0 + gradleft*gradleft + grady*grady);
    143:         coeffright = 1.0/PetscSqrtScalar(1.0 + gradright*gradright + grady*grady);

    145:         f[j][i] = (coeffup*gradup - coeffdown*graddown)*hx + (coeffright*gradright - coeffleft*gradleft)*hy;
    146:
    147:       }

    149:     }
    150:   }
    151:   DAVecRestoreArray((DA)dmmg->dm,localT,&t);
    152:   DAVecRestoreArray((DA)dmmg->dm,F,&f);
    153:   DARestoreLocalVector((DA)dmmg->dm,&localT);
    154:   return(0);
    155: }

    157: PetscErrorCode FormFunctionLocal(DALocalInfo *info,PetscScalar **t,PetscScalar **f,void *ptr)
    158: {
    159:   PetscInt     i,j;
    160:   PetscScalar  hx,hy;
    161:   PetscScalar  gradup,graddown,gradleft,gradright,gradx,grady;
    162:   PetscScalar  coeffup,coeffdown,coeffleft,coeffright;

    165:   hx    = 1.0/(PetscReal)(info->mx-1);  hy    = 1.0/(PetscReal)(info->my-1);
    166:
    167:   /* Evaluate function */
    168:   for (j=info->ys; j<info->ys+info->ym; j++) {
    169:     for (i=info->xs; i<info->xs+info->xm; i++) {

    171:       if (i == 0 || i == info->mx-1 || j == 0 || j == info->my-1) {

    173:         f[j][i] = t[j][i] - (1.0 - (2.0*hx*(PetscReal)i - 1.0)*(2.0*hx*(PetscReal)i - 1.0));
    174:
    175:       } else {

    177:         gradup     = (t[j+1][i] - t[j][i])/hy;
    178:         graddown   = (t[j][i] - t[j-1][i])/hy;
    179:         gradright  = (t[j][i+1] - t[j][i])/hx;
    180:         gradleft   = (t[j][i] - t[j][i-1])/hx;

    182:         gradx      = .5*(t[j][i+1] - t[j][i-1])/hx;
    183:         grady      = .5*(t[j+1][i] - t[j-1][i])/hy;

    185:         coeffup    = 1.0/PetscSqrtScalar(1.0 + gradup*gradup + gradx*gradx);
    186:         coeffdown  = 1.0/PetscSqrtScalar(1.0 + graddown*graddown + gradx*gradx);

    188:         coeffleft  = 1.0/PetscSqrtScalar(1.0 + gradleft*gradleft + grady*grady);
    189:         coeffright = 1.0/PetscSqrtScalar(1.0 + gradright*gradright + grady*grady);

    191:         f[j][i] = (coeffup*gradup - coeffdown*graddown)*hx + (coeffright*gradright - coeffleft*gradleft)*hy;
    192:
    193:       }

    195:     }
    196:   }
    197:   return(0);
    198: }