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 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 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 */ 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 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 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 c99 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