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