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],
4
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))
7
8 which can be reduced to a first order system,
9
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 */
14
15 #include <petscdm.h>
16 #include <petscdmda.h>
17 #include <petscts.h>
18
19 typedef struct {
20 PetscScalar u, v, th;
21 } Field;
22
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 };
31
FormRHSFunction(TS ts,PetscReal t,Vec U,Vec F,PetscCtx ctx)32 static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec U, Vec F, PetscCtx 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;
42
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];
53
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 }
63
FormIFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,PetscCtx ctx)64 static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, PetscCtx 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;
73
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 }
111
112 /* IJacobian - Compute IJacobian = dF/dU + a dF/dUdot */
FormIJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat J,Mat Jpre,PetscCtx ctx)113 PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat J, Mat Jpre, PetscCtx 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;
122
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 };
141
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 };
152
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 };
163
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 }
178
FormInitialSolution(TS ts,Vec U,PetscCtx ctx)179 PetscErrorCode FormInitialSolution(TS ts, Vec U, PetscCtx 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;
189
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 }
206
main(int argc,char ** argv)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;
217
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));
225
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();
240
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));
248
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));
259
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));
265
266 PetscCall(MatDestroy(&J));
267 PetscCall(VecDestroy(&X));
268 PetscCall(TSDestroy(&ts));
269 PetscCall(DMDestroy(&dm));
270 PetscCall(PetscFinalize());
271 return 0;
272 }
273
274 /*TEST
275
276 build:
277 requires: !single !complex
278
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
281
282 TEST*/
283