Actual source code: ex3fastslowsplit.c

      1: static char help[] = "A fast-slow system for testing ARKIMEX.\n";

      3: /*F

      5: \begin{eqnarray}
      6:                  ys' = -2.0*\frac{-1.0+ys^2.0-\cos(t)}{2.0*ys}+0.05*\frac{-2.0+yf^2-\cos(5.0*t)}{2.0*yf}-\frac{\sin(t)}{2.0*ys}\\
      7:                  yf' = 0.05*\frac{-1.0+ys^2-\cos(t)}{2.0*ys}-\frac{-2.0+yf^2-\cos(5.0*t)}{2.0*yf}-5.0*\frac{\sin(5.0*t)}{2.0*yf}\\
      8: \end{eqnarray}

     10: F*/

     12: /*
     13:   This example demonstrates how to use ARKIMEX for solving a fast-slow system. The system is partitioned additively and component-wise at the same time.
     14:   ys stands for the slow component and yf stands for the fast component. On the RHS for yf, only the term -\frac{-2.0+yf^2-\cos(5.0*t)}{2.0*yf} is treated implicitly while the rest is treated explicitly.
     15: */

     17: #include <petscts.h>

     19: typedef struct {
     20:   PetscReal Tf, dt;
     21: } AppCtx;

     23: static PetscErrorCode RHSFunctionslow(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx)
     24: {
     25:   const PetscScalar *u;
     26:   PetscScalar       *f;

     28:   PetscFunctionBegin;
     29:   PetscCall(VecGetArrayRead(U, &u));
     30:   PetscCall(VecGetArray(F, &f));
     31:   f[0] = -2.0 * (-1.0 + u[0] * u[0] - PetscCosScalar(t)) / (2.0 * u[0]) + 0.05 * (-2.0 + u[1] * u[1] - PetscCosScalar(5.0 * t)) / (2.0 * u[1]) - PetscSinScalar(t) / (2.0 * u[0]);
     32:   PetscCall(VecRestoreArrayRead(U, &u));
     33:   PetscCall(VecRestoreArray(F, &f));
     34:   PetscFunctionReturn(PETSC_SUCCESS);
     35: }

     37: static PetscErrorCode RHSFunctionfast(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx)
     38: {
     39:   const PetscScalar *u;
     40:   PetscScalar       *f;

     42:   PetscFunctionBegin;
     43:   PetscCall(VecGetArrayRead(U, &u));
     44:   PetscCall(VecGetArray(F, &f));
     45:   f[0] = 0.05 * (-1.0 + u[0] * u[0] - PetscCosScalar(t)) / (2.0 * u[0]) - 5.0 * PetscSinScalar(5.0 * t) / (2.0 * u[1]);
     46:   PetscCall(VecRestoreArrayRead(U, &u));
     47:   PetscCall(VecRestoreArray(F, &f));
     48:   PetscFunctionReturn(PETSC_SUCCESS);
     49: }

     51: static PetscErrorCode IFunctionfast(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, AppCtx *ctx)
     52: {
     53:   const PetscScalar *u, *udot;
     54:   PetscScalar       *f;

     56:   PetscFunctionBegin;
     57:   PetscCall(VecGetArrayRead(U, &u));
     58:   PetscCall(VecGetArrayRead(Udot, &udot));
     59:   PetscCall(VecGetArray(F, &f));
     60:   f[0] = udot[0] + (-2.0 + u[1] * u[1] - PetscCosScalar(5.0 * t)) / (2.0 * u[1]);
     61:   PetscCall(VecRestoreArrayRead(Udot, &udot));
     62:   PetscCall(VecRestoreArrayRead(U, &u));
     63:   PetscCall(VecRestoreArray(F, &f));
     64:   PetscFunctionReturn(PETSC_SUCCESS);
     65: }

     67: static PetscErrorCode IJacobianfast(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, AppCtx *ctx)
     68: {
     69:   PetscInt           rowcol[] = {0};
     70:   const PetscScalar *u;
     71:   PetscScalar        J[1][1];

     73:   PetscFunctionBeginUser;
     74:   PetscCall(VecGetArrayRead(U, &u));
     75:   J[0][0] = a + (2.0 + PetscCosScalar(5.0 * t)) / (2.0 * u[1] * u[1]) + 0.5;
     76:   PetscCall(MatSetValues(B, 1, rowcol, 1, rowcol, &J[0][0], INSERT_VALUES));
     77:   PetscCall(VecRestoreArrayRead(U, &u));

     79:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
     80:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
     81:   if (A != B) {
     82:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
     83:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
     84:   }
     85:   PetscFunctionReturn(PETSC_SUCCESS);
     86: }

     88: /*
     89:   Define the analytic solution for check method easily
     90: */
     91: static PetscErrorCode sol_true(PetscReal t, Vec U)
     92: {
     93:   PetscScalar *u;

     95:   PetscFunctionBegin;
     96:   PetscCall(VecGetArray(U, &u));
     97:   u[0] = PetscSqrtScalar(1.0 + PetscCosScalar(t));
     98:   u[1] = PetscSqrtScalar(2.0 + PetscCosScalar(5.0 * t));
     99:   PetscCall(VecRestoreArray(U, &u));
    100:   PetscFunctionReturn(PETSC_SUCCESS);
    101: }

    103: int main(int argc, char **argv)
    104: {
    105:   TS           ts; /* ODE integrator */
    106:   Vec          U;  /* solution will be stored here */
    107:   Vec          Utrue;
    108:   Mat          A;
    109:   PetscMPIInt  size;
    110:   AppCtx       ctx;
    111:   PetscScalar *u;
    112:   IS           iss;
    113:   IS           isf;
    114:   PetscInt    *indicess;
    115:   PetscInt    *indicesf;
    116:   PetscInt     n = 2;
    117:   PetscScalar  error, tt;

    119:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    120:      Initialize program
    121:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    122:   PetscFunctionBeginUser;
    123:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    124:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
    125:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");

    127:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    128:     Create index for slow part and fast part
    129:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    130:   PetscCall(PetscMalloc1(1, &indicess));
    131:   indicess[0] = 0;
    132:   PetscCall(PetscMalloc1(1, &indicesf));
    133:   indicesf[0] = 1;
    134:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 1, indicess, PETSC_COPY_VALUES, &iss));
    135:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, 1, indicesf, PETSC_COPY_VALUES, &isf));

    137:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    138:     Create necessary vector
    139:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    140:   PetscCall(VecCreate(PETSC_COMM_WORLD, &U));
    141:   PetscCall(VecSetSizes(U, n, PETSC_DETERMINE));
    142:   PetscCall(VecSetFromOptions(U));
    143:   PetscCall(VecDuplicate(U, &Utrue));
    144:   PetscCall(VecCopy(U, Utrue));

    146:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    147:     Set runtime options
    148:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    149:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ODE options", "");
    150:   ctx.Tf = 0.3;
    151:   ctx.dt = 0.01;
    152:   PetscCall(PetscOptionsScalar("-Tf", "", "", ctx.Tf, &ctx.Tf, NULL));
    153:   PetscCall(PetscOptionsScalar("-dt", "", "", ctx.dt, &ctx.dt, NULL));
    154:   PetscOptionsEnd();

    156:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    157:      Initialize the solution
    158:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    159:   PetscCall(VecGetArray(U, &u));
    160:   u[0] = PetscSqrtScalar(2.0);
    161:   u[1] = PetscSqrtScalar(3.0);
    162:   PetscCall(VecRestoreArray(U, &u));

    164:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    165:      Create timestepping solver context
    166:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    167:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    168:   PetscCall(TSSetType(ts, TSARKIMEX));
    169:   PetscCall(TSARKIMEXSetFastSlowSplit(ts, PETSC_TRUE));

    171:   PetscCall(TSRHSSplitSetIS(ts, "slow", iss));
    172:   PetscCall(TSRHSSplitSetIS(ts, "fast", isf));
    173:   PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, (TSRHSFunctionFn *)RHSFunctionslow, &ctx));
    174:   PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, (TSRHSFunctionFn *)RHSFunctionfast, &ctx));
    175:   PetscCall(TSRHSSplitSetIFunction(ts, "fast", NULL, (TSIFunctionFn *)IFunctionfast, &ctx));
    176:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
    177:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 1, 1));
    178:   PetscCall(MatSetFromOptions(A));
    179:   PetscCall(MatSetUp(A));
    180:   PetscCall(TSRHSSplitSetIJacobian(ts, "fast", A, A, (TSIJacobianFn *)IJacobianfast, &ctx));

    182:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    183:      Set initial conditions
    184:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    185:   PetscCall(TSSetSolution(ts, U));

    187:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    188:      Set solver options
    189:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    190:   PetscCall(TSSetMaxTime(ts, ctx.Tf));
    191:   PetscCall(TSSetTimeStep(ts, ctx.dt));
    192:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
    193:   PetscCall(TSSetFromOptions(ts));

    195:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    196:      Solve linear system
    197:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    198:   PetscCall(TSSolve(ts, U));
    199:   PetscCall(VecView(U, PETSC_VIEWER_STDOUT_WORLD));

    201:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    202:      Check the error of the PETSc solution
    203:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    204:   PetscCall(TSGetTime(ts, &tt));
    205:   PetscCall(sol_true(tt, Utrue));
    206:   PetscCall(VecAXPY(Utrue, -1.0, U));
    207:   PetscCall(VecNorm(Utrue, NORM_2, &error));

    209:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    210:      Print norm2 error
    211:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    212:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "l2 error norm: %g\n", (double)PetscAbsScalar(error)));

    214:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    215:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
    216:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    217:   PetscCall(VecDestroy(&U));
    218:   PetscCall(TSDestroy(&ts));
    219:   PetscCall(VecDestroy(&Utrue));
    220:   PetscCall(ISDestroy(&iss));
    221:   PetscCall(ISDestroy(&isf));
    222:   PetscCall(PetscFree(indicess));
    223:   PetscCall(PetscFree(indicesf));
    224:   PetscCall(MatDestroy(&A));
    225:   PetscCall(PetscFinalize());
    226:   return 0;
    227: }

    229: /*TEST
    230:     build:
    231:       requires: !complex

    233:     test:

    235:     test:
    236:       suffix: 2
    237:       args: -ts_arkimex_initial_guess_extrapolate 1
    238:       output_file: output/ex3fastslowsplit_1.out

    240: TEST*/