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