"Fossies" - the Fresh Open Source Software Archive

Member "petsc-3.12.3/src/ts/examples/tutorials/autodiff/adr_ex1.cxx" (29 Mar 2019, 11284 Bytes) of package /linux/misc/petsc-3.12.3.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file.

    1 static char help[] = "Demonstrates automatic Jacobian generation using ADOL-C for a nonlinear reaction problem from chemistry.\n";
    2 
    3 /*
    4    Concepts: TS^time-dependent nonlinear problems
    5    Concepts: Automatic differentiation using ADOL-C
    6    Processors: 1
    7 */
    8 /*
    9    REQUIRES configuration of PETSc with option --download-adolc.
   10 
   11    For documentation on ADOL-C, see
   12      $PETSC_ARCH/externalpackages/ADOL-C-2.6.0/ADOL-C/doc/adolc-manual.pdf
   13 */
   14 /* ------------------------------------------------------------------------
   15   See ../advection-diffusion-reaction/ex1 for a description of the problem
   16   ------------------------------------------------------------------------- */
   17 #include <petscts.h>
   18 #include "adolc-utils/drivers.cxx"
   19 #include <adolc/adolc.h>
   20 
   21 typedef struct {
   22   PetscScalar k;
   23   Vec         initialsolution;
   24   AdolcCtx    *adctx; /* Automatic differentiation support */
   25 } AppCtx;
   26 
   27 PetscErrorCode IFunctionView(AppCtx *ctx,PetscViewer v)
   28 {
   29   PetscErrorCode ierr;
   30 
   31   PetscFunctionBegin;
   32   ierr = PetscViewerBinaryWrite(v,&ctx->k,1,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr);
   33   PetscFunctionReturn(0);
   34 }
   35 
   36 PetscErrorCode IFunctionLoad(AppCtx **ctx,PetscViewer v)
   37 {
   38   PetscErrorCode ierr;
   39 
   40   PetscFunctionBegin;
   41   ierr = PetscNew(ctx);CHKERRQ(ierr);
   42   ierr = PetscViewerBinaryRead(v,&(*ctx)->k,1,NULL,PETSC_SCALAR);CHKERRQ(ierr);
   43   PetscFunctionReturn(0);
   44 }
   45 
   46 /*
   47   Defines the ODE passed to the ODE solver
   48 */
   49 PetscErrorCode IFunctionPassive(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,AppCtx *ctx)
   50 {
   51   PetscErrorCode    ierr;
   52   PetscScalar       *f;
   53   const PetscScalar *u,*udot;
   54 
   55   PetscFunctionBegin;
   56   /*  The next three lines allow us to access the entries of the vectors directly */
   57   ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr);
   58   ierr = VecGetArrayRead(Udot,&udot);CHKERRQ(ierr);
   59   ierr = VecGetArray(F,&f);CHKERRQ(ierr);
   60   f[0] = udot[0] + ctx->k*u[0]*u[1];
   61   f[1] = udot[1] + ctx->k*u[0]*u[1];
   62   f[2] = udot[2] - ctx->k*u[0]*u[1];
   63   ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
   64   ierr = VecRestoreArrayRead(Udot,&udot);CHKERRQ(ierr);
   65   ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr);
   66   PetscFunctionReturn(0);
   67 }
   68 
   69 /*
   70   'Active' ADOL-C annotated version, marking dependence upon u.
   71 */
   72 PetscErrorCode IFunctionActive1(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,AppCtx *ctx)
   73 {
   74   PetscErrorCode    ierr;
   75   PetscScalar       *f;
   76   const PetscScalar *u,*udot;
   77 
   78   adouble           f_a[3]; /* 'active' double for dependent variables */
   79   adouble           u_a[3]; /* 'active' double for independent variables */
   80 
   81   PetscFunctionBegin;
   82   /*  The next three lines allow us to access the entries of the vectors directly */
   83   ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr);
   84   ierr = VecGetArrayRead(Udot,&udot);CHKERRQ(ierr);
   85   ierr = VecGetArray(F,&f);CHKERRQ(ierr);
   86 
   87   /* Start of active section */
   88   trace_on(1);
   89   u_a[0] <<= u[0]; u_a[1] <<= u[1]; u_a[2] <<= u[2]; /* Mark independence */
   90   f_a[0] = udot[0] + ctx->k*u_a[0]*u_a[1];
   91   f_a[1] = udot[1] + ctx->k*u_a[0]*u_a[1];
   92   f_a[2] = udot[2] - ctx->k*u_a[0]*u_a[1];
   93   f_a[0] >>= f[0]; f_a[1] >>= f[1]; f_a[2] >>= f[2]; /* Mark dependence */
   94   trace_off();
   95   /* End of active section */
   96 
   97   ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
   98   ierr = VecRestoreArrayRead(Udot,&udot);CHKERRQ(ierr);
   99   ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr);
  100   PetscFunctionReturn(0);
  101 }
  102 
  103 /*
  104   'Active' ADOL-C annotated version, marking dependence upon udot.
  105 */
  106 PetscErrorCode IFunctionActive2(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,AppCtx *ctx)
  107 {
  108   PetscErrorCode    ierr;
  109   PetscScalar       *f;
  110   const PetscScalar *u,*udot;
  111 
  112   adouble           f_a[3];    /* 'active' double for dependent variables */
  113   adouble           udot_a[3]; /* 'active' double for independent variables */
  114 
  115   PetscFunctionBegin;
  116   /*  The next three lines allow us to access the entries of the vectors directly */
  117   ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr);
  118   ierr = VecGetArrayRead(Udot,&udot);CHKERRQ(ierr);
  119   ierr = VecGetArray(F,&f);CHKERRQ(ierr);
  120 
  121   /* Start of active section */
  122   trace_on(2);
  123   udot_a[0] <<= udot[0]; udot_a[1] <<= udot[1]; udot_a[2] <<= udot[2]; /* Mark independence */
  124   f_a[0] = udot_a[0] + ctx->k*u[0]*u[1];
  125   f_a[1] = udot_a[1] + ctx->k*u[0]*u[1];
  126   f_a[2] = udot_a[2] - ctx->k*u[0]*u[1];
  127   f_a[0] >>= f[0]; f_a[1] >>= f[1]; f_a[2] >>= f[2];                   /* Mark dependence */
  128   trace_off();
  129   /* End of active section */
  130 
  131   ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
  132   ierr = VecRestoreArrayRead(Udot,&udot);CHKERRQ(ierr);
  133   ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr);
  134   PetscFunctionReturn(0);
  135 }
  136 
  137 /*
  138  Defines the Jacobian of the ODE passed to the ODE solver, using the PETSc-ADOL-C driver for
  139  implicit TS.
  140 */
  141 PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,AppCtx *ctx)
  142 {
  143   PetscErrorCode ierr;
  144   AppCtx         *appctx = (AppCtx*)ctx;
  145   PetscScalar    *u;
  146 
  147   PetscFunctionBegin;
  148   ierr = VecGetArray(U,&u);CHKERRQ(ierr);
  149   ierr = PetscAdolcComputeIJacobian(1,2,A,u,a,appctx->adctx);CHKERRQ(ierr);
  150   ierr = VecRestoreArray(U,&u);CHKERRQ(ierr);
  151   PetscFunctionReturn(0);
  152 }
  153 
  154 /*
  155      Defines the exact (analytic) solution to the ODE
  156 */
  157 static PetscErrorCode Solution(TS ts,PetscReal t,Vec U,AppCtx *ctx)
  158 {
  159   PetscErrorCode    ierr;
  160   const PetscScalar *uinit;
  161   PetscScalar       *u,d0,q;
  162 
  163   PetscFunctionBegin;
  164   ierr = VecGetArrayRead(ctx->initialsolution,&uinit);CHKERRQ(ierr);
  165   ierr = VecGetArray(U,&u);CHKERRQ(ierr);
  166   d0   = uinit[0] - uinit[1];
  167   if (d0 == 0.0) q = ctx->k*t;
  168   else q = (1.0 - PetscExpScalar(-ctx->k*t*d0))/d0;
  169   u[0] = uinit[0]/(1.0 + uinit[1]*q);
  170   u[1] = u[0] - d0;
  171   u[2] = uinit[1] + uinit[2] - u[1];
  172   ierr = VecRestoreArray(U,&u);CHKERRQ(ierr);
  173   ierr = VecRestoreArrayRead(ctx->initialsolution,&uinit);CHKERRQ(ierr);
  174   PetscFunctionReturn(0);
  175 }
  176 
  177 int main(int argc,char **argv)
  178 {
  179   TS             ts;            /* ODE integrator */
  180   Vec            U,Udot,R;      /* solution, derivative, residual */
  181   Mat            A;             /* Jacobian matrix */
  182   PetscErrorCode ierr;
  183   PetscMPIInt    size;
  184   PetscInt       n = 3;
  185   AppCtx         ctx;
  186   AdolcCtx       *adctx;
  187   PetscScalar    *u;
  188   const char     * const names[] = {"U1","U2","U3",NULL};
  189 
  190   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  191      Initialize program
  192      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  193   ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
  194   ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  195   if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");
  196   ierr = PetscNew(&adctx);CHKERRQ(ierr);
  197   adctx->m = n;adctx->n = n;adctx->p = n;
  198   ctx.adctx = adctx;
  199 
  200   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  201     Create necessary matrix and vectors
  202     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  203   ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  204   ierr = MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
  205   ierr = MatSetFromOptions(A);CHKERRQ(ierr);
  206   ierr = MatSetUp(A);CHKERRQ(ierr);
  207 
  208   ierr = MatCreateVecs(A,&U,NULL);CHKERRQ(ierr);
  209 
  210   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  211     Set runtime options
  212     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  213   ctx.k = .9;
  214   ierr  = PetscOptionsGetScalar(NULL,NULL,"-k",&ctx.k,NULL);CHKERRQ(ierr);
  215   ierr  = VecDuplicate(U,&ctx.initialsolution);CHKERRQ(ierr);
  216   ierr  = VecGetArray(ctx.initialsolution,&u);CHKERRQ(ierr);
  217   u[0]  = 1;
  218   u[1]  = .7;
  219   u[2]  = 0;
  220   ierr  = VecRestoreArray(ctx.initialsolution,&u);CHKERRQ(ierr);
  221   ierr  = PetscOptionsGetVec(NULL,NULL,"-initial",ctx.initialsolution,NULL);CHKERRQ(ierr);
  222 
  223   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  224      Create timestepping solver context
  225      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  226   ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
  227   ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);
  228   ierr = TSSetType(ts,TSROSW);CHKERRQ(ierr);
  229   ierr = TSSetIFunction(ts,NULL,(TSIFunction) IFunctionPassive,&ctx);CHKERRQ(ierr);
  230 
  231   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  232      Set initial conditions
  233    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  234   ierr = Solution(ts,0,U,&ctx);CHKERRQ(ierr);
  235   ierr = TSSetSolution(ts,U);CHKERRQ(ierr);
  236 
  237   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  238      Trace just once for each tape
  239      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  240   ierr = VecDuplicate(U,&Udot);CHKERRQ(ierr);
  241   ierr = VecDuplicate(U,&R);CHKERRQ(ierr);
  242   ierr = IFunctionActive1(ts,0.,U,Udot,R,&ctx);CHKERRQ(ierr);
  243   ierr = IFunctionActive2(ts,0.,U,Udot,R,&ctx);CHKERRQ(ierr);
  244   ierr = VecDestroy(&R);CHKERRQ(ierr);
  245   ierr = VecDestroy(&Udot);CHKERRQ(ierr);
  246 
  247   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  248      Set Jacobian
  249      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  250   ierr = TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,&ctx);CHKERRQ(ierr);
  251   ierr = TSSetSolutionFunction(ts,(TSSolutionFunction)Solution,&ctx);CHKERRQ(ierr);
  252 
  253   {
  254     DM   dm;
  255     void *ptr;
  256     ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
  257     ierr = PetscDLSym(NULL,"IFunctionView",&ptr);CHKERRQ(ierr);
  258     ierr = PetscDLSym(NULL,"IFunctionLoad",&ptr);CHKERRQ(ierr);
  259     ierr = DMTSSetIFunctionSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad);CHKERRQ(ierr);
  260     ierr = DMTSSetIJacobianSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad);CHKERRQ(ierr);
  261   }
  262 
  263   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  264      Set solver options
  265    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  266   ierr = TSSetTimeStep(ts,.001);CHKERRQ(ierr);
  267   ierr = TSSetMaxSteps(ts,1000);CHKERRQ(ierr);
  268   ierr = TSSetMaxTime(ts,20.0);CHKERRQ(ierr);
  269   ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);
  270   ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
  271   ierr = TSMonitorLGSetVariableNames(ts,names);CHKERRQ(ierr);
  272 
  273   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  274      Solve nonlinear system
  275      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  276   ierr = TSSolve(ts,U);CHKERRQ(ierr);
  277 
  278   ierr = TSView(ts,PETSC_VIEWER_BINARY_WORLD);CHKERRQ(ierr);
  279 
  280   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  281      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
  282    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  283   ierr = VecDestroy(&ctx.initialsolution);CHKERRQ(ierr);
  284   ierr = MatDestroy(&A);CHKERRQ(ierr);
  285   ierr = VecDestroy(&U);CHKERRQ(ierr);
  286   ierr = TSDestroy(&ts);CHKERRQ(ierr);
  287   ierr = PetscFree(adctx);CHKERRQ(ierr);
  288   ierr = PetscFinalize();
  289   return ierr;
  290 }
  291 
  292 /*TEST
  293 
  294   build:
  295     requires: double !complex adolc
  296 
  297   test:
  298     suffix: 1
  299     args: -ts_max_steps 10 -ts_monitor -ts_adjoint_monitor
  300     output_file: output/adr_ex1_1.out
  301 
  302   test:
  303     suffix: 2
  304     args: -ts_max_steps 1 -snes_test_jacobian
  305     output_file: output/adr_ex1_2.out
  306 
  307 TEST*/