Actual source code: ex34.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  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:   FormInitialSolution(ts, X, &user);
259:   TSSetSolution(ts, X);
260:   VecGetSize(X, &mx);
261:   hx   = 20.0/(PetscReal)(mx-1);
262:   dt   = 0.4 * PetscSqr(hx) / PetscSqr(user.c); /* Diffusive stability limit */
263:   TSSetInitialTimeStep(ts, 0.0, dt);
264:   TSSetFromOptions(ts);

266:   TSSolve(ts, X);
267:   TSGetSolveTime(ts, &ftime);
268:   TSGetTimeStepNumber(ts, &steps);
269:   TSGetConvergedReason(ts, &reason);
270:   PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %D steps\n", TSConvergedReasons[reason], (double)ftime, steps);

272:   MatDestroy(&J);
273:   VecDestroy(&X);
274:   TSDestroy(&ts);
275:   DMDestroy(&dm);
276:   PetscFinalize();
277:   return 0;
278: }