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