Actual source code: spectraladjointassimilation.c

      1: static char help[] = "Solves a simple data assimilation problem with one dimensional advection diffusion equation using TSAdjoint\n\n";

      3: /*

      5:     Not yet tested in parallel

      7: */

      9: /* ------------------------------------------------------------------------

     11:    This program uses the one-dimensional advection-diffusion equation),
     12:        u_t = mu*u_xx - a u_x,
     13:    on the domain 0 <= x <= 1, with periodic boundary conditions

     15:    to demonstrate solving a data assimilation problem of finding the initial conditions
     16:    to produce a given solution at a fixed time.

     18:    The operators are discretized with the spectral element method

     20:   ------------------------------------------------------------------------- */

     22: /*
     23:    Include "petscts.h" so that we can use TS solvers.  Note that this file
     24:    automatically includes:
     25:      petscsys.h       - base PETSc routines   petscvec.h  - vectors
     26:      petscmat.h  - matrices
     27:      petscis.h     - index sets            petscksp.h  - Krylov subspace methods
     28:      petscviewer.h - viewers               petscpc.h   - preconditioners
     29:      petscksp.h   - linear solvers        petscsnes.h - nonlinear solvers
     30: */

     32: #include <petsctao.h>
     33: #include <petscts.h>
     34: #include <petscdt.h>
     35: #include <petscdraw.h>
     36: #include <petscdmda.h>

     38: /*
     39:    User-defined application context - contains data needed by the
     40:    application-provided call-back routines.
     41: */

     43: typedef struct {
     44:   PetscInt   n;       /* number of nodes */
     45:   PetscReal *nodes;   /* GLL nodes */
     46:   PetscReal *weights; /* GLL weights */
     47: } PetscGLL;

     49: typedef struct {
     50:   PetscInt  N;               /* grid points per elements*/
     51:   PetscInt  E;               /* number of elements */
     52:   PetscReal tol_L2, tol_max; /* error norms */
     53:   PetscInt  steps;           /* number of timesteps */
     54:   PetscReal Tend;            /* endtime */
     55:   PetscReal mu;              /* viscosity */
     56:   PetscReal a;               /* advection speed */
     57:   PetscReal L;               /* total length of domain */
     58:   PetscReal Le;
     59:   PetscReal Tadj;
     60: } PetscParam;

     62: typedef struct {
     63:   Vec reference; /* desired end state */
     64:   Vec grid;      /* total grid */
     65:   Vec grad;
     66:   Vec ic;
     67:   Vec curr_sol;
     68:   Vec joe;
     69:   Vec true_solution; /* actual initial conditions for the final solution */
     70: } PetscData;

     72: typedef struct {
     73:   Vec      grid;  /* total grid */
     74:   Vec      mass;  /* mass matrix for total integration */
     75:   Mat      stiff; /* stifness matrix */
     76:   Mat      advec;
     77:   Mat      keptstiff;
     78:   PetscGLL gll;
     79: } PetscSEMOperators;

     81: typedef struct {
     82:   DM                da; /* distributed array data structure */
     83:   PetscSEMOperators SEMop;
     84:   PetscParam        param;
     85:   PetscData         dat;
     86:   TS                ts;
     87:   PetscReal         initial_dt;
     88:   PetscReal        *solutioncoefficients;
     89:   PetscInt          ncoeff;
     90: } AppCtx;

     92: /*
     93:    User-defined routines
     94: */
     95: extern PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
     96: extern PetscErrorCode RHSLaplacian(TS, PetscReal, Vec, Mat, Mat, void *);
     97: extern PetscErrorCode RHSAdvection(TS, PetscReal, Vec, Mat, Mat, void *);
     98: extern PetscErrorCode InitialConditions(Vec, AppCtx *);
     99: extern PetscErrorCode ComputeReference(TS, PetscReal, Vec, AppCtx *);
    100: extern PetscErrorCode MonitorError(Tao, void *);
    101: extern PetscErrorCode MonitorDestroy(void **);
    102: extern PetscErrorCode ComputeSolutionCoefficients(AppCtx *);
    103: extern PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *);
    104: extern PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);

    106: int main(int argc, char **argv)
    107: {
    108:   AppCtx       appctx; /* user-defined application context */
    109:   Tao          tao;
    110:   Vec          u; /* approximate solution vector */
    111:   PetscInt     i, xs, xm, ind, j, lenglob;
    112:   PetscReal    x, *wrk_ptr1, *wrk_ptr2;
    113:   MatNullSpace nsp;

    115:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    116:      Initialize program and set problem parameters
    117:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    118:   PetscFunctionBeginUser;
    119:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));

    121:   /*initialize parameters */
    122:   appctx.param.N     = 10;      /* order of the spectral element */
    123:   appctx.param.E     = 8;       /* number of elements */
    124:   appctx.param.L     = 1.0;     /* length of the domain */
    125:   appctx.param.mu    = 0.00001; /* diffusion coefficient */
    126:   appctx.param.a     = 0.0;     /* advection speed */
    127:   appctx.initial_dt  = 1e-4;
    128:   appctx.param.steps = PETSC_INT_MAX;
    129:   appctx.param.Tend  = 0.01;
    130:   appctx.ncoeff      = 2;

    132:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &appctx.param.N, NULL));
    133:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-E", &appctx.param.E, NULL));
    134:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ncoeff", &appctx.ncoeff, NULL));
    135:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-Tend", &appctx.param.Tend, NULL));
    136:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-mu", &appctx.param.mu, NULL));
    137:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-a", &appctx.param.a, NULL));
    138:   appctx.param.Le = appctx.param.L / appctx.param.E;

    140:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    141:      Create GLL data structures
    142:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    143:   PetscCall(PetscMalloc2(appctx.param.N, &appctx.SEMop.gll.nodes, appctx.param.N, &appctx.SEMop.gll.weights));
    144:   PetscCall(PetscDTGaussLobattoLegendreQuadrature(appctx.param.N, PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA, appctx.SEMop.gll.nodes, appctx.SEMop.gll.weights));
    145:   appctx.SEMop.gll.n = appctx.param.N;
    146:   lenglob            = appctx.param.E * (appctx.param.N - 1);

    148:   /*
    149:      Create distributed array (DMDA) to manage parallel grid and vectors
    150:      and to set up the ghost point communication pattern.  There are E*(Nl-1)+1
    151:      total grid values spread equally among all the processors, except first and last
    152:   */

    154:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, lenglob, 1, 1, NULL, &appctx.da));
    155:   PetscCall(DMSetFromOptions(appctx.da));
    156:   PetscCall(DMSetUp(appctx.da));

    158:   /*
    159:      Extract global and local vectors from DMDA; we use these to store the
    160:      approximate solution.  Then duplicate these for remaining vectors that
    161:      have the same types.
    162:   */

    164:   PetscCall(DMCreateGlobalVector(appctx.da, &u));
    165:   PetscCall(VecDuplicate(u, &appctx.dat.ic));
    166:   PetscCall(VecDuplicate(u, &appctx.dat.true_solution));
    167:   PetscCall(VecDuplicate(u, &appctx.dat.reference));
    168:   PetscCall(VecDuplicate(u, &appctx.SEMop.grid));
    169:   PetscCall(VecDuplicate(u, &appctx.SEMop.mass));
    170:   PetscCall(VecDuplicate(u, &appctx.dat.curr_sol));
    171:   PetscCall(VecDuplicate(u, &appctx.dat.joe));

    173:   PetscCall(DMDAGetCorners(appctx.da, &xs, NULL, NULL, &xm, NULL, NULL));
    174:   PetscCall(DMDAVecGetArray(appctx.da, appctx.SEMop.grid, &wrk_ptr1));
    175:   PetscCall(DMDAVecGetArray(appctx.da, appctx.SEMop.mass, &wrk_ptr2));

    177:   /* Compute function over the locally owned part of the grid */

    179:   xs = xs / (appctx.param.N - 1);
    180:   xm = xm / (appctx.param.N - 1);

    182:   /*
    183:      Build total grid and mass over entire mesh (multi-elemental)
    184:   */

    186:   for (i = xs; i < xs + xm; i++) {
    187:     for (j = 0; j < appctx.param.N - 1; j++) {
    188:       x             = (appctx.param.Le / 2.0) * (appctx.SEMop.gll.nodes[j] + 1.0) + appctx.param.Le * i;
    189:       ind           = i * (appctx.param.N - 1) + j;
    190:       wrk_ptr1[ind] = x;
    191:       wrk_ptr2[ind] = .5 * appctx.param.Le * appctx.SEMop.gll.weights[j];
    192:       if (j == 0) wrk_ptr2[ind] += .5 * appctx.param.Le * appctx.SEMop.gll.weights[j];
    193:     }
    194:   }
    195:   PetscCall(DMDAVecRestoreArray(appctx.da, appctx.SEMop.grid, &wrk_ptr1));
    196:   PetscCall(DMDAVecRestoreArray(appctx.da, appctx.SEMop.mass, &wrk_ptr2));

    198:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    199:    Create matrix data structure; set matrix evaluation routine.
    200:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    201:   PetscCall(DMSetMatrixPreallocateOnly(appctx.da, PETSC_TRUE));
    202:   PetscCall(DMCreateMatrix(appctx.da, &appctx.SEMop.stiff));
    203:   PetscCall(DMCreateMatrix(appctx.da, &appctx.SEMop.advec));

    205:   /*
    206:    For linear problems with a time-dependent f(u,t) in the equation
    207:    u_t = f(u,t), the user provides the discretized right-hand side
    208:    as a time-dependent matrix.
    209:    */
    210:   PetscCall(RHSLaplacian(appctx.ts, 0.0, u, appctx.SEMop.stiff, appctx.SEMop.stiff, &appctx));
    211:   PetscCall(RHSAdvection(appctx.ts, 0.0, u, appctx.SEMop.advec, appctx.SEMop.advec, &appctx));
    212:   PetscCall(MatAXPY(appctx.SEMop.stiff, -1.0, appctx.SEMop.advec, DIFFERENT_NONZERO_PATTERN));
    213:   PetscCall(MatDuplicate(appctx.SEMop.stiff, MAT_COPY_VALUES, &appctx.SEMop.keptstiff));

    215:   /* attach the null space to the matrix, this probably is not needed but does no harm */
    216:   PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, NULL, &nsp));
    217:   PetscCall(MatSetNullSpace(appctx.SEMop.stiff, nsp));
    218:   PetscCall(MatNullSpaceTest(nsp, appctx.SEMop.stiff, NULL));
    219:   PetscCall(MatNullSpaceDestroy(&nsp));

    221:   /* Create the TS solver that solves the ODE and its adjoint; set its options */
    222:   PetscCall(TSCreate(PETSC_COMM_WORLD, &appctx.ts));
    223:   PetscCall(TSSetSolutionFunction(appctx.ts, (PetscErrorCode (*)(TS, PetscReal, Vec, void *))ComputeReference, &appctx));
    224:   PetscCall(TSSetProblemType(appctx.ts, TS_LINEAR));
    225:   PetscCall(TSSetType(appctx.ts, TSRK));
    226:   PetscCall(TSSetDM(appctx.ts, appctx.da));
    227:   PetscCall(TSSetTime(appctx.ts, 0.0));
    228:   PetscCall(TSSetTimeStep(appctx.ts, appctx.initial_dt));
    229:   PetscCall(TSSetMaxSteps(appctx.ts, appctx.param.steps));
    230:   PetscCall(TSSetMaxTime(appctx.ts, appctx.param.Tend));
    231:   PetscCall(TSSetExactFinalTime(appctx.ts, TS_EXACTFINALTIME_MATCHSTEP));
    232:   PetscCall(TSSetTolerances(appctx.ts, 1e-7, NULL, 1e-7, NULL));
    233:   PetscCall(TSSetFromOptions(appctx.ts));
    234:   /* Need to save initial timestep user may have set with -ts_dt so it can be reset for each new TSSolve() */
    235:   PetscCall(TSGetTimeStep(appctx.ts, &appctx.initial_dt));
    236:   PetscCall(TSSetRHSFunction(appctx.ts, NULL, TSComputeRHSFunctionLinear, &appctx));
    237:   PetscCall(TSSetRHSJacobian(appctx.ts, appctx.SEMop.stiff, appctx.SEMop.stiff, TSComputeRHSJacobianConstant, &appctx));
    238:   /*  PetscCall(TSSetRHSFunction(appctx.ts,NULL,RHSFunction,&appctx));
    239:       PetscCall(TSSetRHSJacobian(appctx.ts,appctx.SEMop.stiff,appctx.SEMop.stiff,RHSJacobian,&appctx)); */

    241:   /* Set random initial conditions as initial guess, compute analytic reference solution and analytic (true) initial conditions */
    242:   PetscCall(ComputeSolutionCoefficients(&appctx));
    243:   PetscCall(InitialConditions(appctx.dat.ic, &appctx));
    244:   PetscCall(ComputeReference(appctx.ts, appctx.param.Tend, appctx.dat.reference, &appctx));
    245:   PetscCall(ComputeReference(appctx.ts, 0.0, appctx.dat.true_solution, &appctx));

    247:   /* Set up to save trajectory before TSSetFromOptions() so that TSTrajectory options can be captured */
    248:   PetscCall(TSSetSaveTrajectory(appctx.ts));
    249:   PetscCall(TSSetFromOptions(appctx.ts));

    251:   /* Create TAO solver and set desired solution method  */
    252:   PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
    253:   PetscCall(TaoMonitorSet(tao, MonitorError, &appctx, MonitorDestroy));
    254:   PetscCall(TaoSetType(tao, TAOBQNLS));
    255:   PetscCall(TaoSetSolution(tao, appctx.dat.ic));
    256:   /* Set routine for function and gradient evaluation  */
    257:   PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&appctx));
    258:   /* Check for any TAO command line options  */
    259:   PetscCall(TaoSetTolerances(tao, 1e-8, PETSC_CURRENT, PETSC_CURRENT));
    260:   PetscCall(TaoSetFromOptions(tao));
    261:   PetscCall(TaoSolve(tao));

    263:   PetscCall(TaoDestroy(&tao));
    264:   PetscCall(PetscFree(appctx.solutioncoefficients));
    265:   PetscCall(MatDestroy(&appctx.SEMop.advec));
    266:   PetscCall(MatDestroy(&appctx.SEMop.stiff));
    267:   PetscCall(MatDestroy(&appctx.SEMop.keptstiff));
    268:   PetscCall(VecDestroy(&u));
    269:   PetscCall(VecDestroy(&appctx.dat.ic));
    270:   PetscCall(VecDestroy(&appctx.dat.joe));
    271:   PetscCall(VecDestroy(&appctx.dat.true_solution));
    272:   PetscCall(VecDestroy(&appctx.dat.reference));
    273:   PetscCall(VecDestroy(&appctx.SEMop.grid));
    274:   PetscCall(VecDestroy(&appctx.SEMop.mass));
    275:   PetscCall(VecDestroy(&appctx.dat.curr_sol));
    276:   PetscCall(PetscFree2(appctx.SEMop.gll.nodes, appctx.SEMop.gll.weights));
    277:   PetscCall(DMDestroy(&appctx.da));
    278:   PetscCall(TSDestroy(&appctx.ts));

    280:   /*
    281:      Always call PetscFinalize() before exiting a program.  This routine
    282:        - finalizes the PETSc libraries as well as MPI
    283:        - provides summary and diagnostic information if certain runtime
    284:          options are chosen (e.g., -log_view).
    285:   */
    286:   PetscCall(PetscFinalize());
    287:   return 0;
    288: }

    290: /*
    291:     Computes the coefficients for the analytic solution to the PDE
    292: */
    293: PetscErrorCode ComputeSolutionCoefficients(AppCtx *appctx)
    294: {
    295:   PetscRandom rand;
    296:   PetscInt    i;

    298:   PetscFunctionBegin;
    299:   PetscCall(PetscMalloc1(appctx->ncoeff, &appctx->solutioncoefficients));
    300:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
    301:   PetscCall(PetscRandomSetInterval(rand, .9, 1.0));
    302:   for (i = 0; i < appctx->ncoeff; i++) PetscCall(PetscRandomGetValue(rand, &appctx->solutioncoefficients[i]));
    303:   PetscCall(PetscRandomDestroy(&rand));
    304:   PetscFunctionReturn(PETSC_SUCCESS);
    305: }

    307: /* --------------------------------------------------------------------- */
    308: /*
    309:    InitialConditions - Computes the (random) initial conditions for the Tao optimization solve (these are also initial conditions for the first TSSolve()

    311:    Input Parameter:
    312:    u - uninitialized solution vector (global)
    313:    appctx - user-defined application context

    315:    Output Parameter:
    316:    u - vector with solution at initial time (global)
    317: */
    318: PetscErrorCode InitialConditions(Vec u, AppCtx *appctx)
    319: {
    320:   PetscScalar       *s;
    321:   const PetscScalar *xg;
    322:   PetscInt           i, j, lenglob;
    323:   PetscReal          sum, val;
    324:   PetscRandom        rand;

    326:   PetscFunctionBegin;
    327:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
    328:   PetscCall(PetscRandomSetInterval(rand, .9, 1.0));
    329:   PetscCall(DMDAVecGetArray(appctx->da, u, &s));
    330:   PetscCall(DMDAVecGetArrayRead(appctx->da, appctx->SEMop.grid, (void *)&xg));
    331:   lenglob = appctx->param.E * (appctx->param.N - 1);
    332:   for (i = 0; i < lenglob; i++) {
    333:     s[i] = 0;
    334:     for (j = 0; j < appctx->ncoeff; j++) {
    335:       PetscCall(PetscRandomGetValue(rand, &val));
    336:       s[i] += val * PetscSinScalar(2 * (j + 1) * PETSC_PI * xg[i]);
    337:     }
    338:   }
    339:   PetscCall(PetscRandomDestroy(&rand));
    340:   PetscCall(DMDAVecRestoreArray(appctx->da, u, &s));
    341:   PetscCall(DMDAVecRestoreArrayRead(appctx->da, appctx->SEMop.grid, (void *)&xg));
    342:   /* make sure initial conditions do not contain the constant functions, since with periodic boundary conditions the constant functions introduce a null space */
    343:   PetscCall(VecSum(u, &sum));
    344:   PetscCall(VecShift(u, -sum / lenglob));
    345:   PetscFunctionReturn(PETSC_SUCCESS);
    346: }

    348: /*
    349:    TrueSolution() computes the true solution for the Tao optimization solve which means they are the initial conditions for the objective function.

    351:              InitialConditions() computes the initial conditions for the beginning of the Tao iterations

    353:    Input Parameter:
    354:    u - uninitialized solution vector (global)
    355:    appctx - user-defined application context

    357:    Output Parameter:
    358:    u - vector with solution at initial time (global)
    359: */
    360: PetscErrorCode TrueSolution(Vec u, AppCtx *appctx)
    361: {
    362:   PetscScalar       *s;
    363:   const PetscScalar *xg;
    364:   PetscInt           i, j, lenglob;
    365:   PetscReal          sum;

    367:   PetscFunctionBegin;
    368:   PetscCall(DMDAVecGetArray(appctx->da, u, &s));
    369:   PetscCall(DMDAVecGetArrayRead(appctx->da, appctx->SEMop.grid, (void *)&xg));
    370:   lenglob = appctx->param.E * (appctx->param.N - 1);
    371:   for (i = 0; i < lenglob; i++) {
    372:     s[i] = 0;
    373:     for (j = 0; j < appctx->ncoeff; j++) s[i] += appctx->solutioncoefficients[j] * PetscSinScalar(2 * (j + 1) * PETSC_PI * xg[i]);
    374:   }
    375:   PetscCall(DMDAVecRestoreArray(appctx->da, u, &s));
    376:   PetscCall(DMDAVecRestoreArrayRead(appctx->da, appctx->SEMop.grid, (void *)&xg));
    377:   /* make sure initial conditions do not contain the constant functions, since with periodic boundary conditions the constant functions introduce a null space */
    378:   PetscCall(VecSum(u, &sum));
    379:   PetscCall(VecShift(u, -sum / lenglob));
    380:   PetscFunctionReturn(PETSC_SUCCESS);
    381: }
    382: /* --------------------------------------------------------------------- */
    383: /*
    384:    Sets the desired profile for the final end time

    386:    Input Parameters:
    387:    t - final time
    388:    obj - vector storing the desired profile
    389:    appctx - user-defined application context

    391: */
    392: PetscErrorCode ComputeReference(TS ts, PetscReal t, Vec obj, AppCtx *appctx)
    393: {
    394:   PetscScalar       *s, tc;
    395:   const PetscScalar *xg;
    396:   PetscInt           i, j, lenglob;

    398:   PetscFunctionBegin;
    399:   PetscCall(DMDAVecGetArray(appctx->da, obj, &s));
    400:   PetscCall(DMDAVecGetArrayRead(appctx->da, appctx->SEMop.grid, (void *)&xg));
    401:   lenglob = appctx->param.E * (appctx->param.N - 1);
    402:   for (i = 0; i < lenglob; i++) {
    403:     s[i] = 0;
    404:     for (j = 0; j < appctx->ncoeff; j++) {
    405:       tc = -appctx->param.mu * (j + 1) * (j + 1) * 4.0 * PETSC_PI * PETSC_PI * t;
    406:       s[i] += appctx->solutioncoefficients[j] * PetscSinScalar(2 * (j + 1) * PETSC_PI * (xg[i] + appctx->param.a * t)) * PetscExpReal(tc);
    407:     }
    408:   }
    409:   PetscCall(DMDAVecRestoreArray(appctx->da, obj, &s));
    410:   PetscCall(DMDAVecRestoreArrayRead(appctx->da, appctx->SEMop.grid, (void *)&xg));
    411:   PetscFunctionReturn(PETSC_SUCCESS);
    412: }

    414: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec globalin, Vec globalout, void *ctx)
    415: {
    416:   AppCtx *appctx = (AppCtx *)ctx;

    418:   PetscFunctionBegin;
    419:   PetscCall(MatMult(appctx->SEMop.keptstiff, globalin, globalout));
    420:   PetscFunctionReturn(PETSC_SUCCESS);
    421: }

    423: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec globalin, Mat A, Mat B, void *ctx)
    424: {
    425:   AppCtx *appctx = (AppCtx *)ctx;

    427:   PetscFunctionBegin;
    428:   PetscCall(MatCopy(appctx->SEMop.keptstiff, A, DIFFERENT_NONZERO_PATTERN));
    429:   PetscFunctionReturn(PETSC_SUCCESS);
    430: }

    432: /* --------------------------------------------------------------------- */

    434: /*
    435:    RHSLaplacian -   matrix for diffusion

    437:    Input Parameters:
    438:    ts - the TS context
    439:    t - current time  (ignored)
    440:    X - current solution (ignored)
    441:    dummy - optional user-defined context, as set by TSetRHSJacobian()

    443:    Output Parameters:
    444:    AA - Jacobian matrix
    445:    BB - optionally different matrix from which the preconditioner is built

    447:    Scales by the inverse of the mass matrix (perhaps that should be pulled out)

    449: */
    450: PetscErrorCode RHSLaplacian(TS ts, PetscReal t, Vec X, Mat A, Mat BB, void *ctx)
    451: {
    452:   PetscReal **temp;
    453:   PetscReal   vv;
    454:   AppCtx     *appctx = (AppCtx *)ctx; /* user-defined application context */
    455:   PetscInt    i, xs, xn, l, j;
    456:   PetscInt   *rowsDM;

    458:   PetscFunctionBegin;
    459:   /*
    460:    Creates the element stiffness matrix for the given gll
    461:    */
    462:   PetscCall(PetscGaussLobattoLegendreElementLaplacianCreate(appctx->SEMop.gll.n, appctx->SEMop.gll.nodes, appctx->SEMop.gll.weights, &temp));

    464:   /* scale by the size of the element */
    465:   for (i = 0; i < appctx->param.N; i++) {
    466:     vv = -appctx->param.mu * 2.0 / appctx->param.Le;
    467:     for (j = 0; j < appctx->param.N; j++) temp[i][j] = temp[i][j] * vv;
    468:   }

    470:   PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
    471:   PetscCall(DMDAGetCorners(appctx->da, &xs, NULL, NULL, &xn, NULL, NULL));

    473:   PetscCheck(appctx->param.N - 1 >= 1, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Polynomial order must be at least 2");
    474:   xs = xs / (appctx->param.N - 1);
    475:   xn = xn / (appctx->param.N - 1);

    477:   PetscCall(PetscMalloc1(appctx->param.N, &rowsDM));
    478:   /*
    479:    loop over local elements
    480:    */
    481:   for (j = xs; j < xs + xn; j++) {
    482:     for (l = 0; l < appctx->param.N; l++) rowsDM[l] = 1 + (j - xs) * (appctx->param.N - 1) + l;
    483:     PetscCall(MatSetValuesLocal(A, appctx->param.N, rowsDM, appctx->param.N, rowsDM, &temp[0][0], ADD_VALUES));
    484:   }
    485:   PetscCall(PetscFree(rowsDM));
    486:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    487:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    488:   PetscCall(VecReciprocal(appctx->SEMop.mass));
    489:   PetscCall(MatDiagonalScale(A, appctx->SEMop.mass, 0));
    490:   PetscCall(VecReciprocal(appctx->SEMop.mass));

    492:   PetscCall(PetscGaussLobattoLegendreElementLaplacianDestroy(appctx->SEMop.gll.n, appctx->SEMop.gll.nodes, appctx->SEMop.gll.weights, &temp));
    493:   PetscFunctionReturn(PETSC_SUCCESS);
    494: }

    496: /*
    497:     Almost identical to Laplacian

    499:     Note that the element matrix is NOT scaled by the size of element like the Laplacian term.
    500:  */
    501: PetscErrorCode RHSAdvection(TS ts, PetscReal t, Vec X, Mat A, Mat BB, void *ctx)
    502: {
    503:   PetscReal **temp;
    504:   PetscReal   vv;
    505:   AppCtx     *appctx = (AppCtx *)ctx; /* user-defined application context */
    506:   PetscInt    i, xs, xn, l, j;
    507:   PetscInt   *rowsDM;

    509:   PetscFunctionBegin;
    510:   /*
    511:    Creates the element stiffness matrix for the given gll
    512:    */
    513:   PetscCall(PetscGaussLobattoLegendreElementAdvectionCreate(appctx->SEMop.gll.n, appctx->SEMop.gll.nodes, appctx->SEMop.gll.weights, &temp));

    515:   /* scale by the size of the element */
    516:   for (i = 0; i < appctx->param.N; i++) {
    517:     vv = -appctx->param.a;
    518:     for (j = 0; j < appctx->param.N; j++) temp[i][j] = temp[i][j] * vv;
    519:   }

    521:   PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
    522:   PetscCall(DMDAGetCorners(appctx->da, &xs, NULL, NULL, &xn, NULL, NULL));

    524:   PetscCheck(appctx->param.N - 1 >= 1, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Polynomial order must be at least 2");
    525:   xs = xs / (appctx->param.N - 1);
    526:   xn = xn / (appctx->param.N - 1);

    528:   PetscCall(PetscMalloc1(appctx->param.N, &rowsDM));
    529:   /*
    530:    loop over local elements
    531:    */
    532:   for (j = xs; j < xs + xn; j++) {
    533:     for (l = 0; l < appctx->param.N; l++) rowsDM[l] = 1 + (j - xs) * (appctx->param.N - 1) + l;
    534:     PetscCall(MatSetValuesLocal(A, appctx->param.N, rowsDM, appctx->param.N, rowsDM, &temp[0][0], ADD_VALUES));
    535:   }
    536:   PetscCall(PetscFree(rowsDM));
    537:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    538:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    539:   PetscCall(VecReciprocal(appctx->SEMop.mass));
    540:   PetscCall(MatDiagonalScale(A, appctx->SEMop.mass, 0));
    541:   PetscCall(VecReciprocal(appctx->SEMop.mass));

    543:   PetscCall(PetscGaussLobattoLegendreElementAdvectionDestroy(appctx->SEMop.gll.n, appctx->SEMop.gll.nodes, appctx->SEMop.gll.weights, &temp));
    544:   PetscFunctionReturn(PETSC_SUCCESS);
    545: }

    547: /* ------------------------------------------------------------------ */
    548: /*
    549:    FormFunctionGradient - Evaluates the function and corresponding gradient.

    551:    Input Parameters:
    552:    tao - the Tao context
    553:    ic   - the input vector
    554:    ctx - optional user-defined context, as set when calling TaoSetObjectiveAndGradient()

    556:    Output Parameters:
    557:    f   - the newly evaluated function
    558:    G   - the newly evaluated gradient

    560:    Notes:

    562:           The forward equation is
    563:               M u_t = F(U)
    564:           which is converted to
    565:                 u_t = M^{-1} F(u)
    566:           in the user code since TS has no direct way of providing a mass matrix. The Jacobian of this is
    567:                  M^{-1} J
    568:           where J is the Jacobian of F. Now the adjoint equation is
    569:                 M v_t = J^T v
    570:           but TSAdjoint does not solve this since it can only solve the transposed system for the
    571:           Jacobian the user provided. Hence TSAdjoint solves
    572:                  w_t = J^T M^{-1} w  (where w = M v)
    573:           since there is no way to indicate the mass matrix as a separate entity to TS. Thus one
    574:           must be careful in initializing the "adjoint equation" and using the result. This is
    575:           why
    576:               G = -2 M(u(T) - u_d)
    577:           below (instead of -2(u(T) - u_d)

    579: */
    580: PetscErrorCode FormFunctionGradient(Tao tao, Vec ic, PetscReal *f, Vec G, void *ctx)
    581: {
    582:   AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
    583:   Vec     temp;

    585:   PetscFunctionBegin;
    586:   PetscCall(TSSetTime(appctx->ts, 0.0));
    587:   PetscCall(TSSetStepNumber(appctx->ts, 0));
    588:   PetscCall(TSSetTimeStep(appctx->ts, appctx->initial_dt));
    589:   PetscCall(VecCopy(ic, appctx->dat.curr_sol));

    591:   PetscCall(TSSolve(appctx->ts, appctx->dat.curr_sol));
    592:   PetscCall(VecCopy(appctx->dat.curr_sol, appctx->dat.joe));

    594:   /*     Compute the difference between the current ODE solution and target ODE solution */
    595:   PetscCall(VecWAXPY(G, -1.0, appctx->dat.curr_sol, appctx->dat.reference));

    597:   /*     Compute the objective/cost function   */
    598:   PetscCall(VecDuplicate(G, &temp));
    599:   PetscCall(VecPointwiseMult(temp, G, G));
    600:   PetscCall(VecDot(temp, appctx->SEMop.mass, f));
    601:   PetscCall(VecDestroy(&temp));

    603:   /*     Compute initial conditions for the adjoint integration. See Notes above  */
    604:   PetscCall(VecScale(G, -2.0));
    605:   PetscCall(VecPointwiseMult(G, G, appctx->SEMop.mass));
    606:   PetscCall(TSSetCostGradients(appctx->ts, 1, &G, NULL));

    608:   PetscCall(TSAdjointSolve(appctx->ts));
    609:   /* PetscCall(VecPointwiseDivide(G,G,appctx->SEMop.mass));*/
    610:   PetscFunctionReturn(PETSC_SUCCESS);
    611: }

    613: PetscErrorCode MonitorError(Tao tao, void *ctx)
    614: {
    615:   AppCtx   *appctx = (AppCtx *)ctx;
    616:   Vec       temp, grad;
    617:   PetscReal nrm;
    618:   PetscInt  its;
    619:   PetscReal fct, gnorm;

    621:   PetscFunctionBegin;
    622:   PetscCall(VecDuplicate(appctx->dat.ic, &temp));
    623:   PetscCall(VecWAXPY(temp, -1.0, appctx->dat.ic, appctx->dat.true_solution));
    624:   PetscCall(VecPointwiseMult(temp, temp, temp));
    625:   PetscCall(VecDot(temp, appctx->SEMop.mass, &nrm));
    626:   nrm = PetscSqrtReal(nrm);
    627:   PetscCall(TaoGetGradient(tao, &grad, NULL, NULL));
    628:   PetscCall(VecPointwiseMult(temp, temp, temp));
    629:   PetscCall(VecDot(temp, appctx->SEMop.mass, &gnorm));
    630:   gnorm = PetscSqrtReal(gnorm);
    631:   PetscCall(VecDestroy(&temp));
    632:   PetscCall(TaoGetIterationNumber(tao, &its));
    633:   PetscCall(TaoGetSolutionStatus(tao, NULL, &fct, NULL, NULL, NULL, NULL));
    634:   if (!its) {
    635:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%% Iteration Error Objective Gradient-norm\n"));
    636:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "history = [\n"));
    637:   }
    638:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%3" PetscInt_FMT " %g %g %g\n", its, (double)nrm, (double)fct, (double)gnorm));
    639:   PetscFunctionReturn(PETSC_SUCCESS);
    640: }

    642: PetscErrorCode MonitorDestroy(void **ctx)
    643: {
    644:   PetscFunctionBegin;
    645:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "];\n"));
    646:   PetscFunctionReturn(PETSC_SUCCESS);
    647: }

    649: /*TEST

    651:    build:
    652:      requires: !complex

    654:    test:
    655:      requires: !single
    656:      args: -ts_adapt_dt_max 3.e-3 -E 10 -N 8 -ncoeff 5 -tao_bqnls_mat_lmvm_scale_type none

    658:    test:
    659:      suffix: cn
    660:      requires: !single
    661:      args: -ts_type cn -ts_dt .003 -pc_type lu -E 10 -N 8 -ncoeff 5 -tao_bqnls_mat_lmvm_scale_type none

    663:    test:
    664:      suffix: 2
    665:      requires: !single
    666:      args: -ts_adapt_dt_max 3.e-3 -E 10 -N 8 -ncoeff 5 -a .1 -tao_bqnls_mat_lmvm_scale_type none

    668: TEST*/