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