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;
42: PetscErrorCode ierr;
45: TSGetDM(ts, &dm);
46: DMGetCoordinateDM(dm, &cdm);
47: DMGetCoordinatesLocal(dm, &C);
48: DMDAGetLocalInfo(dm, &info);
49: DMDAVecGetArrayRead(dm, U, (void*)&u);
50: DMDAVecGetArray(dm, F, &f);
51: DMDAVecGetArrayRead(cdm, C, (void*)&x);
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];
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: DMDAVecRestoreArrayRead(dm, U, (void*)&u);
60: DMDAVecRestoreArray(dm, F, &f);
61: DMDAVecRestoreArrayRead(cdm, C, (void*)&x);
62: return(0);
63: }
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;
77: TSGetDM(ts, &dm);
78: DMDAGetLocalInfo(dm, &info);
79: DMGetCoordinateDM(dm, &cdm);
80: DMGetCoordinatesLocal(dm, &C);
81: DMGetLocalVector(dm, &Uloc);
82: DMGlobalToLocalBegin(dm, U, INSERT_VALUES, Uloc);
83: DMGlobalToLocalEnd(dm, U, INSERT_VALUES, Uloc);
84: DMDAVecGetArrayRead(dm, Uloc, &u);
85: DMDAVecGetArrayRead(dm, Udot, &udot);
86: DMDAVecGetArray(dm, F, &f);
87: DMDAVecGetArrayRead(cdm, C, &x);
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: DMDAVecRestoreArrayRead(dm, Uloc, &u);
107: DMDAVecRestoreArrayRead(dm, Udot, &udot);
108: DMDAVecRestoreArray(dm, F, &f);
109: DMDAVecRestoreArrayRead(cdm, C, &x);
110: DMRestoreLocalVector(dm, &Uloc);
111: return(0);
112: }
114: /* IJacobian - Compute IJacobian = dF/dU + a dF/dUdot */
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;
127: TSGetDM(ts, &dm);
128: DMDAGetLocalInfo(dm, &info);
129: DMGetCoordinateDM(dm, &cdm);
130: DMGetCoordinatesLocal(dm, &C);
131: DMDAVecGetArrayRead(dm, U, &u);
132: DMDAVecGetArrayRead(dm, Udot, &udot);
133: DMDAVecGetArrayRead(cdm, C, &x);
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}}};
143: MatSetValuesBlocked(Jpre, 1, &row, 2, col, &vals[0][0][0], INSERT_VALUES);
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}}};
152: MatSetValuesBlocked(Jpre, 1, &row, 2, col, &vals[0][0][0], INSERT_VALUES);
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}}};
161: MatSetValuesBlocked(Jpre, 1, &row, 3, col, &vals[0][0][0], INSERT_VALUES);
162: }
163: }
164: DMDAVecRestoreArrayRead(dm, U, &u);
165: DMDAVecRestoreArrayRead(dm, Udot, &udot);
166: DMDAVecRestoreArrayRead(cdm, C, &x);
167: MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY);
168: MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY);
169: if (J != Jpre) {
170: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
171: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
172: }
173: return(0);
174: }
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;
189: TSGetDM(ts, &dm);
190: DMGetCoordinateDM(dm, &cdm);
191: DMGetCoordinatesLocal(dm, &C);
192: DMDAGetLocalInfo(dm, &info);
193: DMDAVecGetArray(dm, U, &u);
194: DMDAVecGetArrayRead(cdm, C, &x);
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: DMDAVecRestoreArray(dm, U, &u);
201: DMDAVecRestoreArrayRead(cdm, C, &x);
202: return(0);
203: }
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;
217: PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr;
218: DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 11, 3, 1, NULL, &dm);
219: DMSetFromOptions(dm);
220: DMSetUp(dm);
221: DMDASetUniformCoordinates(dm, 0.0, 20.0, 0.0, 0.0, 0.0, 0.0);
222: DMCreateGlobalVector(dm, &X);
224: 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: PetscOptionsReal("-epsilon", "Inverse of seismic ratio", "", user.epsilon, &user.epsilon, NULL);
232: PetscOptionsReal("-gamma", "Wave frequency for interblock coupling", "", user.gamma, &user.gamma, NULL);
233: PetscOptionsReal("-gamma_tilde", "Wave frequency for plate coupling", "", user.gammaTilde, &user.gammaTilde, NULL);
234: PetscOptionsReal("-xi", "Interblock spring constant", "", user.xi, &user.xi, NULL);
235: PetscOptionsReal("-c", "Wavespeed", "", user.c, &user.c, NULL);
236: }
237: PetscOptionsEnd();
239: TSCreate(PETSC_COMM_WORLD, &ts);
240: TSSetDM(ts, dm);
241: TSSetRHSFunction(ts, NULL, FormRHSFunction, &user);
242: TSSetIFunction(ts, NULL, FormIFunction, &user);
243: DMSetMatType(dm, MATAIJ);
244: DMCreateMatrix(dm, &J);
245: TSSetIJacobian(ts, J, J, FormIJacobian, &user);
247: ftime = 800.0;
248: TSSetMaxTime(ts,ftime);
249: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
250: FormInitialSolution(ts, X, &user);
251: TSSetSolution(ts, X);
252: VecGetSize(X, &mx);
253: hx = 20.0/(PetscReal)(mx-1);
254: dt = 0.4 * PetscSqr(hx) / PetscSqr(user.c); /* Diffusive stability limit */
255: TSSetTimeStep(ts,dt);
256: TSSetFromOptions(ts);
258: TSSolve(ts, X);
259: TSGetSolveTime(ts, &ftime);
260: TSGetStepNumber(ts, &steps);
261: TSGetConvergedReason(ts, &reason);
262: PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %D steps\n", TSConvergedReasons[reason], (double)ftime, steps);
264: MatDestroy(&J);
265: VecDestroy(&X);
266: TSDestroy(&ts);
267: DMDestroy(&dm);
268: PetscFinalize();
269: return ierr;
270: }
272: /*TEST
274: build:
275: requires: !single !complex
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
280: TEST*/