Actual source code: nladic.c

      1: #define PETSCMAT_DLL
      2: /*
      3:     ADIC based nonlinear operator object that can be used with FAS

      5:     This does not really belong in the matrix directories but since it
      6:     was cloned off of Mat_DAAD I'm leaving it here until I have a better place

      8: */

     10:  #include petscda.h
     11:  #include petscsys.h
     13: #include "adic/ad_utils.h"

     16:  #include src/dm/da/daimpl.h

     18: struct NLF_DAAD {
     19:   DA         da;
     20:   void       *ctx;
     21:   Vec        residual;
     22:   int        newton_its;
     23: };

     27: /*
     28:       Solves the one dimensional equation using Newton's method
     29: */
     32: PetscErrorCode NLFNewton_DAAD(NLF A,DALocalInfo *info,MatStencil *stencil,void *ad_vu,PetscScalar *ad_vustart,int nI,int gI,PetscScalar residual)
     33: {
     35:   PetscInt       cnt = A->newton_its;
     36:   PetscScalar    ad_f[2],J,f;

     39:   ad_vustart[1+2*gI] = 1.0;

     41:   do {
     42:     /* compute the function and Jacobian */
     43:     (*A->da->adicmf_lfi)(info,stencil,ad_vu,ad_f,A->ctx);
     44:     J    = -ad_f[1];
     45:     f    = -ad_f[0] + residual;
     46:     ad_vustart[2*gI] =  ad_vustart[2*gI] - f/J;
     47:   } while (--cnt > 0 && PetscAbsScalar(f) > 1.e-14);

     49:   ad_vustart[1+2*gI] = 0.0;
     50:   return(0);
     51: }

     54: /*
     55:         Nonlinear relax on all the equations with an initial guess in xx
     56: */
     60: PetscErrorCode PETSCMAT_DLLEXPORT NLFRelax_DAAD(NLF A,MatSORType flag,int its,Vec xx)
     61: {
     63:   PetscInt       j,gtdof,nI,gI;
     64:   PetscScalar    *avu,*av,*ad_vustart,*residual;
     65:   Vec            localxx;
     66:   DALocalInfo    info;
     67:   MatStencil     stencil;
     68:   void*          *ad_vu;

     71:   if (its <= 0) SETERRQ1(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D positive",its);

     73:   DAGetLocalVector(A->da,&localxx);
     74:   /* get space for derivative object.  */
     75:   DAGetAdicMFArray(A->da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,&gtdof);
     76:   VecGetArray(A->residual,&residual);

     79:   /* tell ADIC we will be computing one dimensional Jacobians */
     80:   PetscADResetIndep();
     81:   PetscADIncrementTotalGradSize(1);
     82:   PetscADSetIndepDone();

     84:   DAGetLocalInfo(A->da,&info);
     85:   while (its--) {

     87:     /* get initial solution properly ghosted */
     88:     DAGlobalToLocalBegin(A->da,xx,INSERT_VALUES,localxx);
     89:     DAGlobalToLocalEnd(A->da,xx,INSERT_VALUES,localxx);

     91:     /* copy input vector into derivative object */
     92:     VecGetArray(localxx,&avu);
     93:     for (j=0; j<gtdof; j++) {
     94:       ad_vustart[2*j]   = avu[j];
     95:       ad_vustart[2*j+1] = 0.0;
     96:     }
     97:     VecRestoreArray(localxx,&avu);

     99:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){
    100:       nI = 0;
    101:       for (stencil.k = info.zs; stencil.k<info.zs+info.zm; stencil.k++) {
    102:         for (stencil.j = info.ys; stencil.j<info.ys+info.ym; stencil.j++) {
    103:           for (stencil.i = info.xs; stencil.i<info.xs+info.xm; stencil.i++) {
    104:             for (stencil.c = 0; stencil.c<info.dof; stencil.c++) {
    105:               gI   = stencil.c + (stencil.i - info.gxs)*info.dof + (stencil.j - info.gys)*info.dof*info.gxm + (stencil.k - info.gzs)*info.dof*info.gxm*info.gym;
    106:               NLFNewton_DAAD(A,&info,&stencil,ad_vu,ad_vustart,nI,gI,residual[nI]);
    107:               nI++;
    108:             }
    109:           }
    110:         }
    111:       }
    112:     }
    113:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){
    114:       nI = info.dof*info.xm*info.ym*info.zm - 1;
    115:       for (stencil.k = info.zs+info.zm-1; stencil.k>=info.zs; stencil.k--) {
    116:         for (stencil.j = info.ys+info.ym-1; stencil.j>=info.ys; stencil.j--) {
    117:           for (stencil.i = info.xs+info.xm-1; stencil.i>=info.xs; stencil.i--) {
    118:             for (stencil.c = info.dof-1; stencil.c>=0; stencil.c--) {
    119:               gI   = stencil.c + (stencil.i - info.gxs)*info.dof + (stencil.j - info.gys)*info.dof*info.gxm + (stencil.k - info.gzs)*info.dof*info.gxm*info.gym;
    120:               NLFNewton_DAAD(A,&info,&stencil,ad_vu,ad_vustart,nI,gI,residual[nI]);
    121:               nI--;
    122:             }
    123:           }
    124:         }
    125:       }
    126:     }

    128:     /* copy solution back into ghosted vector from derivative object */
    129:     VecGetArray(localxx,&av);
    130:     for (j=0; j<gtdof; j++) {
    131:       av[j] = ad_vustart[2*j];
    132:     }
    133:     VecRestoreArray(localxx,&av);
    134:     /* stick relaxed solution back into global solution */
    135:     DALocalToGlobal(A->da,localxx,INSERT_VALUES,xx);
    136:   }

    139:   VecRestoreArray(A->residual,&residual);
    140:   DARestoreLocalVector(A->da,&localxx);
    141:   DARestoreAdicMFArray(A->da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,&gtdof);
    142:   return(0);
    143: }

    149: PetscErrorCode NLFDestroy_DAAD(NLF A)
    150: {

    154:   DADestroy(A->da);
    155:   PetscFree(A);
    156:   return(0);
    157: }

    163: PetscErrorCode PETSCMAT_DLLEXPORT NLFDAADSetDA_DAAD(NLF A,DA da)
    164: {

    168:   A->da = da;
    169:   PetscObjectReference((PetscObject)da);
    170:   return(0);
    171: }

    177: PetscErrorCode PETSCMAT_DLLEXPORT NLFDAADSetNewtonIterations_DAAD(NLF A,int its)
    178: {
    180:   A->newton_its = its;
    181:   return(0);
    182: }

    188: PetscErrorCode PETSCMAT_DLLEXPORT NLFDAADSetResidual_DAAD(NLF A,Vec residual)
    189: {
    191:   A->residual = residual;
    192:   return(0);
    193: }

    200: PetscErrorCode PETSCMAT_DLLEXPORT NLFDAADSetCtx_DAAD(NLF A,void *ctx)
    201: {
    203:   A->ctx = ctx;
    204:   return(0);
    205: }

    211: PetscErrorCode PETSCMAT_DLLEXPORT NLFCreate_DAAD(NLF *A)
    212: {

    216:   PetscNew(struct NLF_DAAD,A);
    217:   (*A)->da         = 0;
    218:   (*A)->ctx        = 0;
    219:   (*A)->newton_its = 2;
    220:   return(0);
    221: }