Actual source code: ex34.c

      1: static const char help[] = "An elastic wave equation driven by Dieterich-Ruina friction\n";
      2: /*
      3: This whole derivation comes from Erickson, Birnir, and Lavallee [2010]. The model comes from the continuum limit in Carlson and Langer [1989],

      5:   u_{tt}   = c^2 u_{xx} - \tilde\gamma^2 u - (\gamma^2 / \xi) (\theta + \ln(u_t + 1))
      6:   \theta_t = -(u_t + 1) (\theta + (1 + \epsilon) \ln(u_t +1))

      8: which can be reduced to a first order system,

     10:   u_t      = v
     11:   v_t      = c^2 u_{xx} - \tilde\gamma^2 u - (\gamma^2 / \xi)(\theta + ln(v + 1)))
     12:   \theta_t = -(v + 1) (\theta + (1 + \epsilon) \ln(v+1))
     13: */

     15: #include <petscdm.h>
     16: #include <petscdmda.h>
     17: #include <petscts.h>

     19: typedef struct {
     20:   PetscScalar u, v, th;
     21: } Field;

     23: typedef struct _User *User;
     24: struct _User {
     25:   PetscReal epsilon;    /* inverse of seismic ratio, B-A / A */
     26:   PetscReal gamma;      /* wave frequency for interblock coupling */
     27:   PetscReal gammaTilde; /* wave frequency for coupling to plate */
     28:   PetscReal xi;         /* interblock spring constant */
     29:   PetscReal c;          /* wavespeed */
     30: };

     32: static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
     33: {
     34:   User               user = (User)ctx;
     35:   DM                 dm, cdm;
     36:   DMDALocalInfo      info;
     37:   Vec                C;
     38:   Field             *f;
     39:   const Field       *u;
     40:   const PetscScalar *x;
     41:   PetscInt           i;

     43:   PetscFunctionBeginUser;
     44:   PetscCall(TSGetDM(ts, &dm));
     45:   PetscCall(DMGetCoordinateDM(dm, &cdm));
     46:   PetscCall(DMGetCoordinatesLocal(dm, &C));
     47:   PetscCall(DMDAGetLocalInfo(dm, &info));
     48:   PetscCall(DMDAVecGetArrayRead(dm, U, (void *)&u));
     49:   PetscCall(DMDAVecGetArray(dm, F, &f));
     50:   PetscCall(DMDAVecGetArrayRead(cdm, C, (void *)&x));
     51:   for (i = info.xs; i < info.xs + info.xm; ++i) {
     52:     const PetscScalar hx = i + 1 == info.xs + info.xm ? x[i] - x[i - 1] : x[i + 1] - x[i];

     54:     f[i].u  = hx * (u[i].v);
     55:     f[i].v  = -hx * (PetscSqr(user->gammaTilde) * u[i].u + (PetscSqr(user->gamma) / user->xi) * (u[i].th + PetscLogScalar(u[i].v + 1)));
     56:     f[i].th = -hx * (u[i].v + 1) * (u[i].th + (1 + user->epsilon) * PetscLogScalar(u[i].v + 1));
     57:   }
     58:   PetscCall(DMDAVecRestoreArrayRead(dm, U, (void *)&u));
     59:   PetscCall(DMDAVecRestoreArray(dm, F, &f));
     60:   PetscCall(DMDAVecRestoreArrayRead(cdm, C, (void *)&x));
     61:   PetscFunctionReturn(PETSC_SUCCESS);
     62: }

     64: static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
     65: {
     66:   User          user = (User)ctx;
     67:   DM            dm, cdm;
     68:   DMDALocalInfo info;
     69:   Vec           Uloc, C;
     70:   Field        *u, *udot, *f;
     71:   PetscScalar  *x;
     72:   PetscInt      i;

     74:   PetscFunctionBeginUser;
     75:   PetscCall(TSGetDM(ts, &dm));
     76:   PetscCall(DMDAGetLocalInfo(dm, &info));
     77:   PetscCall(DMGetCoordinateDM(dm, &cdm));
     78:   PetscCall(DMGetCoordinatesLocal(dm, &C));
     79:   PetscCall(DMGetLocalVector(dm, &Uloc));
     80:   PetscCall(DMGlobalToLocalBegin(dm, U, INSERT_VALUES, Uloc));
     81:   PetscCall(DMGlobalToLocalEnd(dm, U, INSERT_VALUES, Uloc));
     82:   PetscCall(DMDAVecGetArrayRead(dm, Uloc, &u));
     83:   PetscCall(DMDAVecGetArrayRead(dm, Udot, &udot));
     84:   PetscCall(DMDAVecGetArray(dm, F, &f));
     85:   PetscCall(DMDAVecGetArrayRead(cdm, C, &x));
     86:   for (i = info.xs; i < info.xs + info.xm; ++i) {
     87:     if (i == 0) {
     88:       const PetscScalar hx = x[i + 1] - x[i];
     89:       f[i].u               = hx * udot[i].u;
     90:       f[i].v               = hx * udot[i].v - PetscSqr(user->c) * (u[i + 1].u - u[i].u) / hx;
     91:       f[i].th              = hx * udot[i].th;
     92:     } else if (i == info.mx - 1) {
     93:       const PetscScalar hx = x[i] - x[i - 1];
     94:       f[i].u               = hx * udot[i].u;
     95:       f[i].v               = hx * udot[i].v - PetscSqr(user->c) * (u[i - 1].u - u[i].u) / hx;
     96:       f[i].th              = hx * udot[i].th;
     97:     } else {
     98:       const PetscScalar hx = x[i + 1] - x[i];
     99:       f[i].u               = hx * udot[i].u;
    100:       f[i].v               = hx * udot[i].v - PetscSqr(user->c) * (u[i - 1].u - 2. * u[i].u + u[i + 1].u) / hx;
    101:       f[i].th              = hx * udot[i].th;
    102:     }
    103:   }
    104:   PetscCall(DMDAVecRestoreArrayRead(dm, Uloc, &u));
    105:   PetscCall(DMDAVecRestoreArrayRead(dm, Udot, &udot));
    106:   PetscCall(DMDAVecRestoreArray(dm, F, &f));
    107:   PetscCall(DMDAVecRestoreArrayRead(cdm, C, &x));
    108:   PetscCall(DMRestoreLocalVector(dm, &Uloc));
    109:   PetscFunctionReturn(PETSC_SUCCESS);
    110: }

    112: /* IJacobian - Compute IJacobian = dF/dU + a dF/dUdot */
    113: PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat J, Mat Jpre, void *ctx)
    114: {
    115:   User          user = (User)ctx;
    116:   DM            dm, cdm;
    117:   DMDALocalInfo info;
    118:   Vec           C;
    119:   Field        *u, *udot;
    120:   PetscScalar  *x;
    121:   PetscInt      i;

    123:   PetscFunctionBeginUser;
    124:   PetscCall(TSGetDM(ts, &dm));
    125:   PetscCall(DMDAGetLocalInfo(dm, &info));
    126:   PetscCall(DMGetCoordinateDM(dm, &cdm));
    127:   PetscCall(DMGetCoordinatesLocal(dm, &C));
    128:   PetscCall(DMDAVecGetArrayRead(dm, U, &u));
    129:   PetscCall(DMDAVecGetArrayRead(dm, Udot, &udot));
    130:   PetscCall(DMDAVecGetArrayRead(cdm, C, &x));
    131:   for (i = info.xs; i < info.xs + info.xm; ++i) {
    132:     if (i == 0) {
    133:       const PetscScalar hx  = x[i + 1] - x[i];
    134:       const PetscInt    row = i, col[] = {i, i + 1};
    135:       const PetscScalar dxx0 = PetscSqr(user->c) / hx, dxxR = -PetscSqr(user->c) / hx;
    136:       const PetscScalar vals[3][2][3] = {
    137:         {{a * hx, 0, 0},        {0, 0, 0}   },
    138:         {{0, a * hx + dxx0, 0}, {0, dxxR, 0}},
    139:         {{0, 0, a * hx},        {0, 0, 0}   }
    140:       };

    142:       PetscCall(MatSetValuesBlocked(Jpre, 1, &row, 2, col, &vals[0][0][0], INSERT_VALUES));
    143:     } else if (i == info.mx - 1) {
    144:       const PetscScalar hx  = x[i + 1] - x[i];
    145:       const PetscInt    row = i, col[] = {i - 1, i};
    146:       const PetscScalar dxxL = -PetscSqr(user->c) / hx, dxx0 = PetscSqr(user->c) / hx;
    147:       const PetscScalar vals[3][2][3] = {
    148:         {{0, 0, 0},    {a * hx, 0, 0}       },
    149:         {{0, dxxL, 0}, {0, a * hx + dxx0, 0}},
    150:         {{0, 0, 0},    {0, 0, a * hx}       }
    151:       };

    153:       PetscCall(MatSetValuesBlocked(Jpre, 1, &row, 2, col, &vals[0][0][0], INSERT_VALUES));
    154:     } else {
    155:       const PetscScalar hx  = x[i + 1] - x[i];
    156:       const PetscInt    row = i, col[] = {i - 1, i, i + 1};
    157:       const PetscScalar dxxL = -PetscSqr(user->c) / hx, dxx0 = 2. * PetscSqr(user->c) / hx, dxxR = -PetscSqr(user->c) / hx;
    158:       const PetscScalar vals[3][3][3] = {
    159:         {{0, 0, 0},    {a * hx, 0, 0},        {0, 0, 0}   },
    160:         {{0, dxxL, 0}, {0, a * hx + dxx0, 0}, {0, dxxR, 0}},
    161:         {{0, 0, 0},    {0, 0, a * hx},        {0, 0, 0}   }
    162:       };

    164:       PetscCall(MatSetValuesBlocked(Jpre, 1, &row, 3, col, &vals[0][0][0], INSERT_VALUES));
    165:     }
    166:   }
    167:   PetscCall(DMDAVecRestoreArrayRead(dm, U, &u));
    168:   PetscCall(DMDAVecRestoreArrayRead(dm, Udot, &udot));
    169:   PetscCall(DMDAVecRestoreArrayRead(cdm, C, &x));
    170:   PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
    171:   PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
    172:   if (J != Jpre) {
    173:     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
    174:     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
    175:   }
    176:   PetscFunctionReturn(PETSC_SUCCESS);
    177: }

    179: PetscErrorCode FormInitialSolution(TS ts, Vec U, void *ctx)
    180: {
    181:   /* User            user = (User) ctx; */
    182:   DM              dm, cdm;
    183:   DMDALocalInfo   info;
    184:   Vec             C;
    185:   Field          *u;
    186:   PetscScalar    *x;
    187:   const PetscReal sigma = 1.0;
    188:   PetscInt        i;

    190:   PetscFunctionBeginUser;
    191:   PetscCall(TSGetDM(ts, &dm));
    192:   PetscCall(DMGetCoordinateDM(dm, &cdm));
    193:   PetscCall(DMGetCoordinatesLocal(dm, &C));
    194:   PetscCall(DMDAGetLocalInfo(dm, &info));
    195:   PetscCall(DMDAVecGetArray(dm, U, &u));
    196:   PetscCall(DMDAVecGetArrayRead(cdm, C, &x));
    197:   for (i = info.xs; i < info.xs + info.xm; ++i) {
    198:     u[i].u  = 1.5 * PetscExpScalar(-PetscSqr(x[i] - 10) / PetscSqr(sigma));
    199:     u[i].v  = 0.0;
    200:     u[i].th = 0.0;
    201:   }
    202:   PetscCall(DMDAVecRestoreArray(dm, U, &u));
    203:   PetscCall(DMDAVecRestoreArrayRead(cdm, C, &x));
    204:   PetscFunctionReturn(PETSC_SUCCESS);
    205: }

    207: int main(int argc, char **argv)
    208: {
    209:   DM                dm;
    210:   TS                ts;
    211:   Vec               X;
    212:   Mat               J;
    213:   PetscInt          steps, mx;
    214:   PetscReal         ftime, hx, dt;
    215:   TSConvergedReason reason;
    216:   struct _User      user;

    218:   PetscFunctionBeginUser;
    219:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    220:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 11, 3, 1, NULL, &dm));
    221:   PetscCall(DMSetFromOptions(dm));
    222:   PetscCall(DMSetUp(dm));
    223:   PetscCall(DMDASetUniformCoordinates(dm, 0.0, 20.0, 0.0, 0.0, 0.0, 0.0));
    224:   PetscCall(DMCreateGlobalVector(dm, &X));

    226:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Dynamic Friction Options", "");
    227:   {
    228:     user.epsilon    = 0.1;
    229:     user.gamma      = 0.5;
    230:     user.gammaTilde = 0.5;
    231:     user.xi         = 0.5;
    232:     user.c          = 0.5;
    233:     PetscCall(PetscOptionsReal("-epsilon", "Inverse of seismic ratio", "", user.epsilon, &user.epsilon, NULL));
    234:     PetscCall(PetscOptionsReal("-gamma", "Wave frequency for interblock coupling", "", user.gamma, &user.gamma, NULL));
    235:     PetscCall(PetscOptionsReal("-gamma_tilde", "Wave frequency for plate coupling", "", user.gammaTilde, &user.gammaTilde, NULL));
    236:     PetscCall(PetscOptionsReal("-xi", "Interblock spring constant", "", user.xi, &user.xi, NULL));
    237:     PetscCall(PetscOptionsReal("-c", "Wavespeed", "", user.c, &user.c, NULL));
    238:   }
    239:   PetscOptionsEnd();

    241:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    242:   PetscCall(TSSetDM(ts, dm));
    243:   PetscCall(TSSetRHSFunction(ts, NULL, FormRHSFunction, &user));
    244:   PetscCall(TSSetIFunction(ts, NULL, FormIFunction, &user));
    245:   PetscCall(DMSetMatType(dm, MATAIJ));
    246:   PetscCall(DMCreateMatrix(dm, &J));
    247:   PetscCall(TSSetIJacobian(ts, J, J, FormIJacobian, &user));

    249:   ftime = 800.0;
    250:   PetscCall(TSSetMaxTime(ts, ftime));
    251:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
    252:   PetscCall(FormInitialSolution(ts, X, &user));
    253:   PetscCall(TSSetSolution(ts, X));
    254:   PetscCall(VecGetSize(X, &mx));
    255:   hx = 20.0 / (PetscReal)(mx - 1);
    256:   dt = 0.4 * PetscSqr(hx) / PetscSqr(user.c); /* Diffusive stability limit */
    257:   PetscCall(TSSetTimeStep(ts, dt));
    258:   PetscCall(TSSetFromOptions(ts));

    260:   PetscCall(TSSolve(ts, X));
    261:   PetscCall(TSGetSolveTime(ts, &ftime));
    262:   PetscCall(TSGetStepNumber(ts, &steps));
    263:   PetscCall(TSGetConvergedReason(ts, &reason));
    264:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, steps));

    266:   PetscCall(MatDestroy(&J));
    267:   PetscCall(VecDestroy(&X));
    268:   PetscCall(TSDestroy(&ts));
    269:   PetscCall(DMDestroy(&dm));
    270:   PetscCall(PetscFinalize());
    271:   return 0;
    272: }

    274: /*TEST

    276:     build:
    277:       requires: !single !complex

    279:     test:
    280:       TODO: broken, was not nightly tested, SNES solve eventually fails, -snes_test_jacobian indicates Jacobian is wrong, but even -snes_mf_operator fails

    282: TEST*/