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,void * ctx)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;
42 PetscErrorCode ierr;
43
44 PetscFunctionBeginUser;
45 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr);
46 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);
47 ierr = DMGetCoordinatesLocal(dm, &C);CHKERRQ(ierr);
48 ierr = DMDAGetLocalInfo(dm, &info);CHKERRQ(ierr);
49 ierr = DMDAVecGetArrayRead(dm, U, (void*)&u);CHKERRQ(ierr);
50 ierr = DMDAVecGetArray(dm, F, &f);CHKERRQ(ierr);
51 ierr = DMDAVecGetArrayRead(cdm, C, (void*)&x);CHKERRQ(ierr);
52 for (i = info.xs; i < info.xs+info.xm; ++i) {
53 const PetscScalar hx = i+1 == info.xs+info.xm ? x[i] - x[i-1] : x[i+1] - x[i];
54
55 f[i].u = hx*(u[i].v);
56 f[i].v = -hx*(PetscSqr(user->gammaTilde)*u[i].u + (PetscSqr(user->gamma) / user->xi)*(u[i].th + PetscLogScalar(u[i].v + 1)));
57 f[i].th = -hx*(u[i].v + 1)*(u[i].th + (1 + user->epsilon)*PetscLogScalar(u[i].v + 1));
58 }
59 ierr = DMDAVecRestoreArrayRead(dm, U, (void*)&u);CHKERRQ(ierr);
60 ierr = DMDAVecRestoreArray(dm, F, &f);CHKERRQ(ierr);
61 ierr = DMDAVecRestoreArrayRead(cdm, C, (void*)&x);CHKERRQ(ierr);
62 PetscFunctionReturn(0);
63 }
64
FormIFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void * ctx)65 static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
66 {
67 User user = (User) ctx;
68 DM dm, cdm;
69 DMDALocalInfo info;
70 Vec Uloc, C;
71 Field *u, *udot, *f;
72 PetscScalar *x;
73 PetscInt i;
74 PetscErrorCode ierr;
75
76 PetscFunctionBeginUser;
77 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr);
78 ierr = DMDAGetLocalInfo(dm, &info);CHKERRQ(ierr);
79 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);
80 ierr = DMGetCoordinatesLocal(dm, &C);CHKERRQ(ierr);
81 ierr = DMGetLocalVector(dm, &Uloc);CHKERRQ(ierr);
82 ierr = DMGlobalToLocalBegin(dm, U, INSERT_VALUES, Uloc);CHKERRQ(ierr);
83 ierr = DMGlobalToLocalEnd(dm, U, INSERT_VALUES, Uloc);CHKERRQ(ierr);
84 ierr = DMDAVecGetArrayRead(dm, Uloc, &u);CHKERRQ(ierr);
85 ierr = DMDAVecGetArrayRead(dm, Udot, &udot);CHKERRQ(ierr);
86 ierr = DMDAVecGetArray(dm, F, &f);CHKERRQ(ierr);
87 ierr = DMDAVecGetArrayRead(cdm, C, &x);CHKERRQ(ierr);
88 for (i = info.xs; i < info.xs+info.xm; ++i) {
89 if (i == 0) {
90 const PetscScalar hx = x[i+1] - x[i];
91 f[i].u = hx * udot[i].u;
92 f[i].v = hx * udot[i].v - PetscSqr(user->c) * (u[i+1].u - u[i].u) / hx;
93 f[i].th = hx * udot[i].th;
94 } else if (i == info.mx-1) {
95 const PetscScalar hx = x[i] - x[i-1];
96 f[i].u = hx * udot[i].u;
97 f[i].v = hx * udot[i].v - PetscSqr(user->c) * (u[i-1].u - u[i].u) / hx;
98 f[i].th = hx * udot[i].th;
99 } else {
100 const PetscScalar hx = x[i+1] - x[i];
101 f[i].u = hx * udot[i].u;
102 f[i].v = hx * udot[i].v - PetscSqr(user->c) * (u[i-1].u - 2.*u[i].u + u[i+1].u) / hx;
103 f[i].th = hx * udot[i].th;
104 }
105 }
106 ierr = DMDAVecRestoreArrayRead(dm, Uloc, &u);CHKERRQ(ierr);
107 ierr = DMDAVecRestoreArrayRead(dm, Udot, &udot);CHKERRQ(ierr);
108 ierr = DMDAVecRestoreArray(dm, F, &f);CHKERRQ(ierr);
109 ierr = DMDAVecRestoreArrayRead(cdm, C, &x);CHKERRQ(ierr);
110 ierr = DMRestoreLocalVector(dm, &Uloc);CHKERRQ(ierr);
111 PetscFunctionReturn(0);
112 }
113
114 /* IJacobian - Compute IJacobian = dF/dU + a dF/dUdot */
FormIJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat J,Mat Jpre,void * ctx)115 PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat J, Mat Jpre, void *ctx)
116 {
117 User user = (User) ctx;
118 DM dm, cdm;
119 DMDALocalInfo info;
120 Vec C;
121 Field *u, *udot;
122 PetscScalar *x;
123 PetscInt i;
124 PetscErrorCode ierr;
125
126 PetscFunctionBeginUser;
127 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr);
128 ierr = DMDAGetLocalInfo(dm, &info);CHKERRQ(ierr);
129 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);
130 ierr = DMGetCoordinatesLocal(dm, &C);CHKERRQ(ierr);
131 ierr = DMDAVecGetArrayRead(dm, U, &u);CHKERRQ(ierr);
132 ierr = DMDAVecGetArrayRead(dm, Udot, &udot);CHKERRQ(ierr);
133 ierr = DMDAVecGetArrayRead(cdm, C, &x);CHKERRQ(ierr);
134 for (i = info.xs; i < info.xs+info.xm; ++i) {
135 if (i == 0) {
136 const PetscScalar hx = x[i+1] - x[i];
137 const PetscInt row = i, col[] = {i,i+1};
138 const PetscScalar dxx0 = PetscSqr(user->c)/hx,dxxR = -PetscSqr(user->c)/hx;
139 const PetscScalar vals[3][2][3] = {{{a*hx, 0,0},{0,0, 0}},
140 {{0,a*hx+dxx0,0},{0,dxxR,0}},
141 {{0,0, a*hx},{0,0, 0}}};
142
143 ierr = MatSetValuesBlocked(Jpre, 1, &row, 2, col, &vals[0][0][0], INSERT_VALUES);CHKERRQ(ierr);
144 } else if (i == info.mx-1) {
145 const PetscScalar hx = x[i+1] - x[i];
146 const PetscInt row = i, col[] = {i-1,i};
147 const PetscScalar dxxL = -PetscSqr(user->c)/hx, dxx0 = PetscSqr(user->c)/hx;
148 const PetscScalar vals[3][2][3] = {{{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 ierr = MatSetValuesBlocked(Jpre, 1, &row, 2, col, &vals[0][0][0], INSERT_VALUES);CHKERRQ(ierr);
153 } else {
154 const PetscScalar hx = x[i+1] - x[i];
155 const PetscInt row = i, col[] = {i-1,i,i+1};
156 const PetscScalar dxxL = -PetscSqr(user->c)/hx, dxx0 = 2.*PetscSqr(user->c)/hx,dxxR = -PetscSqr(user->c)/hx;
157 const PetscScalar vals[3][3][3] = {{{0,0, 0},{a*hx, 0,0},{0,0, 0}},
158 {{0,dxxL,0},{0,a*hx+dxx0,0},{0,dxxR,0}},
159 {{0,0, 0},{0,0, a*hx},{0,0, 0}}};
160
161 ierr = MatSetValuesBlocked(Jpre, 1, &row, 3, col, &vals[0][0][0], INSERT_VALUES);CHKERRQ(ierr);
162 }
163 }
164 ierr = DMDAVecRestoreArrayRead(dm, U, &u);CHKERRQ(ierr);
165 ierr = DMDAVecRestoreArrayRead(dm, Udot, &udot);CHKERRQ(ierr);
166 ierr = DMDAVecRestoreArrayRead(cdm, C, &x);CHKERRQ(ierr);
167 ierr = MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
168 ierr = MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
169 if (J != Jpre) {
170 ierr = MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
171 ierr = MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
172 }
173 PetscFunctionReturn(0);
174 }
175
FormInitialSolution(TS ts,Vec U,void * ctx)176 PetscErrorCode FormInitialSolution(TS ts, Vec U, void *ctx)
177 {
178 /* User user = (User) ctx; */
179 DM dm, cdm;
180 DMDALocalInfo info;
181 Vec C;
182 Field *u;
183 PetscScalar *x;
184 const PetscReal sigma = 1.0;
185 PetscInt i;
186 PetscErrorCode ierr;
187
188 PetscFunctionBeginUser;
189 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr);
190 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);
191 ierr = DMGetCoordinatesLocal(dm, &C);CHKERRQ(ierr);
192 ierr = DMDAGetLocalInfo(dm, &info);CHKERRQ(ierr);
193 ierr = DMDAVecGetArray(dm, U, &u);CHKERRQ(ierr);
194 ierr = DMDAVecGetArrayRead(cdm, C, &x);CHKERRQ(ierr);
195 for (i = info.xs; i < info.xs+info.xm; ++i) {
196 u[i].u = 1.5 * PetscExpScalar(-PetscSqr(x[i] - 10)/PetscSqr(sigma));
197 u[i].v = 0.0;
198 u[i].th = 0.0;
199 }
200 ierr = DMDAVecRestoreArray(dm, U, &u);CHKERRQ(ierr);
201 ierr = DMDAVecRestoreArrayRead(cdm, C, &x);CHKERRQ(ierr);
202 PetscFunctionReturn(0);
203 }
204
main(int argc,char ** argv)205 int main(int argc, char **argv)
206 {
207 DM dm;
208 TS ts;
209 Vec X;
210 Mat J;
211 PetscInt steps, mx;
212 PetscReal ftime, hx, dt;
213 TSConvergedReason reason;
214 struct _User user;
215 PetscErrorCode ierr;
216
217 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
218 ierr = DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 11, 3, 1, NULL, &dm);CHKERRQ(ierr);
219 ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
220 ierr = DMSetUp(dm);CHKERRQ(ierr);
221 ierr = DMDASetUniformCoordinates(dm, 0.0, 20.0, 0.0, 0.0, 0.0, 0.0);CHKERRQ(ierr);
222 ierr = DMCreateGlobalVector(dm, &X);CHKERRQ(ierr);
223
224 ierr = PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Dynamic Friction Options", "");
225 {
226 user.epsilon = 0.1;
227 user.gamma = 0.5;
228 user.gammaTilde = 0.5;
229 user.xi = 0.5;
230 user.c = 0.5;
231 ierr = PetscOptionsReal("-epsilon", "Inverse of seismic ratio", "", user.epsilon, &user.epsilon, NULL);CHKERRQ(ierr);
232 ierr = PetscOptionsReal("-gamma", "Wave frequency for interblock coupling", "", user.gamma, &user.gamma, NULL);CHKERRQ(ierr);
233 ierr = PetscOptionsReal("-gamma_tilde", "Wave frequency for plate coupling", "", user.gammaTilde, &user.gammaTilde, NULL);CHKERRQ(ierr);
234 ierr = PetscOptionsReal("-xi", "Interblock spring constant", "", user.xi, &user.xi, NULL);CHKERRQ(ierr);
235 ierr = PetscOptionsReal("-c", "Wavespeed", "", user.c, &user.c, NULL);CHKERRQ(ierr);
236 }
237 ierr = PetscOptionsEnd();CHKERRQ(ierr);
238
239 ierr = TSCreate(PETSC_COMM_WORLD, &ts);CHKERRQ(ierr);
240 ierr = TSSetDM(ts, dm);CHKERRQ(ierr);
241 ierr = TSSetRHSFunction(ts, NULL, FormRHSFunction, &user);CHKERRQ(ierr);
242 ierr = TSSetIFunction(ts, NULL, FormIFunction, &user);CHKERRQ(ierr);
243 ierr = DMSetMatType(dm, MATAIJ);CHKERRQ(ierr);
244 ierr = DMCreateMatrix(dm, &J);CHKERRQ(ierr);
245 ierr = TSSetIJacobian(ts, J, J, FormIJacobian, &user);CHKERRQ(ierr);
246
247 ftime = 800.0;
248 ierr = TSSetMaxTime(ts,ftime);CHKERRQ(ierr);
249 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);
250 ierr = FormInitialSolution(ts, X, &user);CHKERRQ(ierr);
251 ierr = TSSetSolution(ts, X);CHKERRQ(ierr);
252 ierr = VecGetSize(X, &mx);CHKERRQ(ierr);
253 hx = 20.0/(PetscReal)(mx-1);
254 dt = 0.4 * PetscSqr(hx) / PetscSqr(user.c); /* Diffusive stability limit */
255 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr);
256 ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
257
258 ierr = TSSolve(ts, X);CHKERRQ(ierr);
259 ierr = TSGetSolveTime(ts, &ftime);CHKERRQ(ierr);
260 ierr = TSGetStepNumber(ts, &steps);CHKERRQ(ierr);
261 ierr = TSGetConvergedReason(ts, &reason);CHKERRQ(ierr);
262 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %D steps\n", TSConvergedReasons[reason], (double)ftime, steps);CHKERRQ(ierr);
263
264 ierr = MatDestroy(&J);CHKERRQ(ierr);
265 ierr = VecDestroy(&X);CHKERRQ(ierr);
266 ierr = TSDestroy(&ts);CHKERRQ(ierr);
267 ierr = DMDestroy(&dm);CHKERRQ(ierr);
268 ierr = PetscFinalize();
269 return ierr;
270 }
271
272 /*TEST
273
274 build:
275 requires: !single !complex
276
277 test:
278 TODO: broken, was not nightly tested, SNES solve eventually fails, -snes_test_jacobian indicates Jacobian is wrong, but even -snes_mf_operator fails
279
280 TEST*/
281