Actual source code: ex74.c

      1: static char help[] = "Solves the constant-coefficient 1D heat equation \n\
      2: with an Implicit Runge-Kutta method using MatKAIJ.                  \n\
      3:                                                                     \n\
      4:     du      d^2 u                                                   \n\
      5:     --  = a ----- ; 0 <= x <= 1;                                    \n\
      6:     dt      dx^2                                                    \n\
      7:                                                                     \n\
      8:   with periodic boundary conditions                                 \n\
      9:                                                                     \n\
     10: 2nd order central discretization in space:                          \n\
     11:                                                                     \n\
     12:    [ d^2 u ]     u_{i+1} - 2u_i + u_{i-1}                           \n\
     13:    [ ----- ]  =  ------------------------                           \n\
     14:    [ dx^2  ]i              h^2                                      \n\
     15:                                                                     \n\
     16:     i = grid index;    h = x_{i+1}-x_i (Uniform)                    \n\
     17:     0 <= i < n         h = 1.0/n                                    \n\
     18:                                                                     \n\
     19: Thus,                                                               \n\
     20:                                                                     \n\
     21:    du                                                               \n\
     22:    --  = Ju;  J = (a/h^2) tridiagonal(1,-2,1)_n                     \n\
     23:    dt                                                               \n\
     24:                                                                     \n\
     25: This example is a TS version of the KSP ex74.c tutorial.            \n";

     27: #include <petscts.h>

     29: typedef enum {
     30:   PHYSICS_DIFFUSION,
     31:   PHYSICS_ADVECTION
     32: } PhysicsType;
     33: const char *const PhysicsTypes[] = {"DIFFUSION", "ADVECTION", "PhysicsType", "PHYSICS_", NULL};

     35: typedef struct Context {
     36:   PetscReal   a;          /* diffusion coefficient      */
     37:   PetscReal   xmin, xmax; /* domain bounds              */
     38:   PetscInt    imax;       /* number of grid points      */
     39:   PhysicsType physics_type;
     40: } UserContext;

     42: static PetscErrorCode ExactSolution(Vec, void *, PetscReal);
     43: static PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);

     45: int main(int argc, char **argv)
     46: {
     47:   TS          ts;
     48:   Mat         A;
     49:   Vec         u, uex;
     50:   UserContext ctxt;
     51:   PetscReal   err, ftime;

     53:   PetscFunctionBeginUser;
     54:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
     55:   /* default value */
     56:   ctxt.a            = 0.1;
     57:   ctxt.xmin         = 0.0;
     58:   ctxt.xmax         = 1.0;
     59:   ctxt.imax         = 40;
     60:   ctxt.physics_type = PHYSICS_DIFFUSION;

     62:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "IRK options", "");
     63:   PetscCall(PetscOptionsReal("-a", "diffusion coefficient", "<1.0>", ctxt.a, &ctxt.a, NULL));
     64:   PetscCall(PetscOptionsInt("-imax", "grid size", "<20>", ctxt.imax, &ctxt.imax, NULL));
     65:   PetscCall(PetscOptionsReal("-xmin", "xmin", "<0.0>", ctxt.xmin, &ctxt.xmin, NULL));
     66:   PetscCall(PetscOptionsReal("-xmax", "xmax", "<1.0>", ctxt.xmax, &ctxt.xmax, NULL));
     67:   PetscCall(PetscOptionsEnum("-physics_type", "Type of process to discretize", "", PhysicsTypes, (PetscEnum)ctxt.physics_type, (PetscEnum *)&ctxt.physics_type, NULL));
     68:   PetscOptionsEnd();

     70:   /* allocate and initialize solution vector and exact solution */
     71:   PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
     72:   PetscCall(VecSetSizes(u, PETSC_DECIDE, ctxt.imax));
     73:   PetscCall(VecSetFromOptions(u));
     74:   PetscCall(VecDuplicate(u, &uex));
     75:   /* initial solution */
     76:   PetscCall(ExactSolution(u, &ctxt, 0.0));

     78:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
     79:   PetscCall(MatSetType(A, MATAIJ));
     80:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, ctxt.imax, ctxt.imax));
     81:   PetscCall(MatSetUp(A));

     83:   /* Create and set options for TS */
     84:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
     85:   PetscCall(TSSetProblemType(ts, TS_LINEAR));
     86:   PetscCall(TSSetTimeStep(ts, 0.125));
     87:   PetscCall(TSSetSolution(ts, u));
     88:   PetscCall(TSSetMaxSteps(ts, 10));
     89:   PetscCall(TSSetMaxTime(ts, 1.0));
     90:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
     91:   PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, &ctxt));
     92:   PetscCall(RHSJacobian(ts, 0, u, A, A, &ctxt));
     93:   PetscCall(TSSetRHSJacobian(ts, A, A, TSComputeRHSJacobianConstant, &ctxt));
     94:   PetscCall(TSSetFromOptions(ts));
     95:   PetscCall(TSSolve(ts, u));

     97:   PetscCall(TSGetSolveTime(ts, &ftime));
     98:   /* exact   solution */
     99:   PetscCall(ExactSolution(uex, &ctxt, ftime));

    101:   /* Calculate error in final solution */
    102:   PetscCall(VecAYPX(uex, -1.0, u));
    103:   PetscCall(VecNorm(uex, NORM_2, &err));
    104:   err = PetscSqrtReal(err * err / ((PetscReal)ctxt.imax));
    105:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "L2 norm of the numerical error = %g (time=%g)\n", (double)err, (double)ftime));

    107:   /* Free up memory */
    108:   PetscCall(TSDestroy(&ts));
    109:   PetscCall(MatDestroy(&A));
    110:   PetscCall(VecDestroy(&uex));
    111:   PetscCall(VecDestroy(&u));
    112:   PetscCall(PetscFinalize());
    113:   return 0;
    114: }

    116: PetscErrorCode ExactSolution(Vec u, void *c, PetscReal t)
    117: {
    118:   UserContext *ctxt = (UserContext *)c;
    119:   PetscInt     i, is, ie;
    120:   PetscScalar *uarr;
    121:   PetscReal    x, dx, a = ctxt->a, pi = PETSC_PI;

    123:   PetscFunctionBeginUser;
    124:   dx = (ctxt->xmax - ctxt->xmin) / ((PetscReal)ctxt->imax);
    125:   PetscCall(VecGetOwnershipRange(u, &is, &ie));
    126:   PetscCall(VecGetArray(u, &uarr));
    127:   for (i = is; i < ie; i++) {
    128:     x = i * dx;
    129:     switch (ctxt->physics_type) {
    130:     case PHYSICS_DIFFUSION:
    131:       uarr[i - is] = PetscExpScalar(-4.0 * pi * pi * a * t) * PetscSinScalar(2 * pi * x);
    132:       break;
    133:     case PHYSICS_ADVECTION:
    134:       uarr[i - is] = PetscSinScalar(2 * pi * (x - a * t));
    135:       break;
    136:     default:
    137:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for physics type %s", PhysicsTypes[ctxt->physics_type]);
    138:     }
    139:   }
    140:   PetscCall(VecRestoreArray(u, &uarr));
    141:   PetscFunctionReturn(PETSC_SUCCESS);
    142: }

    144: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat Jpre, void *ctx)
    145: {
    146:   UserContext *user = (UserContext *)ctx;
    147:   PetscInt     matis, matie, i;
    148:   PetscReal    dx, dx2;

    150:   PetscFunctionBeginUser;
    151:   dx  = (user->xmax - user->xmin) / ((PetscReal)user->imax);
    152:   dx2 = dx * dx;
    153:   PetscCall(MatGetOwnershipRange(J, &matis, &matie));
    154:   for (i = matis; i < matie; i++) {
    155:     PetscScalar values[3];
    156:     PetscInt    col[3];
    157:     switch (user->physics_type) {
    158:     case PHYSICS_DIFFUSION:
    159:       values[0] = user->a * 1.0 / dx2;
    160:       values[1] = -user->a * 2.0 / dx2;
    161:       values[2] = user->a * 1.0 / dx2;
    162:       break;
    163:     case PHYSICS_ADVECTION:
    164:       values[0] = user->a * .5 / dx;
    165:       values[1] = 0.;
    166:       values[2] = -user->a * .5 / dx;
    167:       break;
    168:     default:
    169:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for physics type %s", PhysicsTypes[user->physics_type]);
    170:     }
    171:     /* periodic boundaries */
    172:     if (i == 0) {
    173:       col[0] = user->imax - 1;
    174:       col[1] = i;
    175:       col[2] = i + 1;
    176:     } else if (i == user->imax - 1) {
    177:       col[0] = i - 1;
    178:       col[1] = i;
    179:       col[2] = 0;
    180:     } else {
    181:       col[0] = i - 1;
    182:       col[1] = i;
    183:       col[2] = i + 1;
    184:     }
    185:     PetscCall(MatSetValues(J, 1, &i, 3, col, values, INSERT_VALUES));
    186:   }
    187:   PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
    188:   PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
    189:   PetscFunctionReturn(PETSC_SUCCESS);
    190: }

    192: /*TEST

    194:   test:
    195:     requires: double
    196:     suffix: 1
    197:     nsize: {{1 2}}
    198:     args: -ts_max_steps 5 -ts_monitor -ksp_monitor_short -pc_type pbjacobi -ksp_atol 1e-6 -ts_type irk -ts_irk_nstages 2

    200:   test:
    201:     requires: double
    202:     suffix: 2
    203:     args: -ts_max_steps 5 -ts_monitor -ksp_monitor_short -pc_type pbjacobi -ksp_atol 1e-6 -ts_type irk -ts_irk_nstages 3

    205:   testset:
    206:     requires: hpddm
    207:     args: -ts_max_steps 5 -ts_monitor -ksp_monitor_short -pc_type pbjacobi -ksp_atol 1e-4 -ts_type irk -ts_irk_nstages 3 -ksp_view_final_residual -ksp_hpddm_type gcrodr -ksp_type hpddm
    208:     test:
    209:       suffix: 3
    210:       requires: double
    211:       args: -ksp_hpddm_precision {{single double}shared output}
    212:     test:
    213:       suffix: 3_single
    214:       requires: single
    215:       args: -ksp_hpddm_precision {{single double}shared output}
    216:     test:
    217:       suffix: 3___float128
    218:       requires: __float128
    219:       output_file: output/ex74_3.out
    220:       args: -ksp_hpddm_precision {{double quadruple}shared output}
    221: TEST*/