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*/