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;

 44:   TSGetDM(ts, &dm);
 45:   DMGetCoordinateDM(dm, &cdm);
 46:   DMGetCoordinatesLocal(dm, &C);
 47:   DMDAGetLocalInfo(dm, &info);
 48:   DMDAVecGetArrayRead(dm, U, (void *)&u);
 49:   DMDAVecGetArray(dm, F, &f);
 50:   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:   DMDAVecRestoreArrayRead(dm, U, (void *)&u);
 59:   DMDAVecRestoreArray(dm, F, &f);
 60:   DMDAVecRestoreArrayRead(cdm, C, (void *)&x);
 61:   return 0;
 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;

 75:   TSGetDM(ts, &dm);
 76:   DMDAGetLocalInfo(dm, &info);
 77:   DMGetCoordinateDM(dm, &cdm);
 78:   DMGetCoordinatesLocal(dm, &C);
 79:   DMGetLocalVector(dm, &Uloc);
 80:   DMGlobalToLocalBegin(dm, U, INSERT_VALUES, Uloc);
 81:   DMGlobalToLocalEnd(dm, U, INSERT_VALUES, Uloc);
 82:   DMDAVecGetArrayRead(dm, Uloc, &u);
 83:   DMDAVecGetArrayRead(dm, Udot, &udot);
 84:   DMDAVecGetArray(dm, F, &f);
 85:   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:   DMDAVecRestoreArrayRead(dm, Uloc, &u);
105:   DMDAVecRestoreArrayRead(dm, Udot, &udot);
106:   DMDAVecRestoreArray(dm, F, &f);
107:   DMDAVecRestoreArrayRead(cdm, C, &x);
108:   DMRestoreLocalVector(dm, &Uloc);
109:   return 0;
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;

124:   TSGetDM(ts, &dm);
125:   DMDAGetLocalInfo(dm, &info);
126:   DMGetCoordinateDM(dm, &cdm);
127:   DMGetCoordinatesLocal(dm, &C);
128:   DMDAVecGetArrayRead(dm, U, &u);
129:   DMDAVecGetArrayRead(dm, Udot, &udot);
130:   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:       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:       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:       MatSetValuesBlocked(Jpre, 1, &row, 3, col, &vals[0][0][0], INSERT_VALUES);
165:     }
166:   }
167:   DMDAVecRestoreArrayRead(dm, U, &u);
168:   DMDAVecRestoreArrayRead(dm, Udot, &udot);
169:   DMDAVecRestoreArrayRead(cdm, C, &x);
170:   MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY);
171:   MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY);
172:   if (J != Jpre) {
173:     MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
174:     MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
175:   }
176:   return 0;
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;

191:   TSGetDM(ts, &dm);
192:   DMGetCoordinateDM(dm, &cdm);
193:   DMGetCoordinatesLocal(dm, &C);
194:   DMDAGetLocalInfo(dm, &info);
195:   DMDAVecGetArray(dm, U, &u);
196:   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:   DMDAVecRestoreArray(dm, U, &u);
203:   DMDAVecRestoreArrayRead(cdm, C, &x);
204:   return 0;
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;

219:   PetscInitialize(&argc, &argv, NULL, help);
220:   DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 11, 3, 1, NULL, &dm);
221:   DMSetFromOptions(dm);
222:   DMSetUp(dm);
223:   DMDASetUniformCoordinates(dm, 0.0, 20.0, 0.0, 0.0, 0.0, 0.0);
224:   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:     PetscOptionsReal("-epsilon", "Inverse of seismic ratio", "", user.epsilon, &user.epsilon, NULL);
234:     PetscOptionsReal("-gamma", "Wave frequency for interblock coupling", "", user.gamma, &user.gamma, NULL);
235:     PetscOptionsReal("-gamma_tilde", "Wave frequency for plate coupling", "", user.gammaTilde, &user.gammaTilde, NULL);
236:     PetscOptionsReal("-xi", "Interblock spring constant", "", user.xi, &user.xi, NULL);
237:     PetscOptionsReal("-c", "Wavespeed", "", user.c, &user.c, NULL);
238:   }
239:   PetscOptionsEnd();

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

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

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

266:   MatDestroy(&J);
267:   VecDestroy(&X);
268:   TSDestroy(&ts);
269:   DMDestroy(&dm);
270:   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*/