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