Lines Matching +full:test +full:- +full:experimental
1 static const char help[] = "Steady-state 2D subduction flow, pressure and temperature solver.\n\
3 ---------------------------------ex30 help---------------------------------\n\
4 -OPTION <DEFAULT> = (UNITS) DESCRIPTION.\n\n\
5 -width <320> = (km) width of domain.\n\
6 -depth <300> = (km) depth of domain.\n\
7 -slab_dip <45> = (degrees) dip angle of the slab (determines the grid aspect ratio).\n\
8 -lid_depth <35> = (km) depth of the static conductive lid.\n\
9 -fault_depth <35> = (km) depth of slab-wedge mechanical coupling\n\
12 -ni <82> = grid cells in x-direction. (nj adjusts to accommodate\n\
13 the slab dip & depth). DO NOT USE -da_grid_x option!!!\n\
14 -ivisc <3> = rheology option.\n\
15 0 --- constant viscosity.\n\
16 1 --- olivine diffusion creep rheology (T&P-dependent, newtonian).\n\
17 2 --- olivine dislocation creep rheology (T&P-dependent, non-newtonian).\n\
18 3 --- Full mantle rheology, combination of 1 & 2.\n\
20 -slab_velocity <5> = (cm/year) convergence rate of slab into subduction zone.\n\
21 -slab_age <50> = (million yrs) age of slab for thermal profile boundary condition.\n\
22 -lid_age <50> = (million yrs) age of lid for thermal profile boundary condition.\n\
25 ---------------------------------ex30 help---------------------------------\n";
27 /*F-----------------------------------------------------------------------
35 -\nabla P + \nabla \cdot [\eta (\nabla v + \nabla v^T)] & = & 0 \\
37 dT/dt + \nabla \cdot (vT) - 1/Pe \triangle^2(T) & = & 0 \\
42 & = & \hbox{diffusion creep (T,P-dependent) } \hbox{if ivisc} ==1 \\
43 & = & \hbox{dislocation creep (T,P,v-dependent)} \hbox{if ivisc} ==2 \\
48 -------$w_{ij}$------
49 $u_{i-1j}$ $P,T_{ij}$ $u_{ij}$
50 ------$w_{ij-1}$-----
52 ------------------------------------------------------------------------F*/
117 /* Post-processing & misc */
126 AppCtx *user; /* user-defined work context */ in main()
135 PetscCall(PetscOptionsSetValue(NULL, "-file", "ex30_output")); in main()
136 PetscCall(PetscOptionsSetValue(NULL, "-snes_monitor_short", NULL)); in main()
137 PetscCall(PetscOptionsSetValue(NULL, "-snes_max_it", "20")); in main()
138 PetscCall(PetscOptionsSetValue(NULL, "-ksp_max_it", "1500")); in main()
139 PetscCall(PetscOptionsSetValue(NULL, "-ksp_gmres_restart", "300")); in main()
144 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - in main()
146 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ in main()
150 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - in main()
151 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ in main()
157 PetscCall(DMDASetFieldName(da, 0, "x-velocity")); in main()
158 PetscCall(DMDASetFieldName(da, 1, "y-velocity")); in main()
162 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - in main()
164 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ in main()
166 user->param = ¶m; in main()
167 user->grid = &grid; in main()
169 PetscCall(DMCreateGlobalVector(da, &user->Xguess)); in main()
171 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - in main()
173 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ in main()
180 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - in main()
182 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ in main()
186 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - in main()
188 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ in main()
191 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - in main()
193 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ in main()
194 PetscCall(VecDestroy(&user->Xguess)); in main()
195 PetscCall(VecDestroy(&user->x)); in main()
214 Parameter *param = user->param; in UpdateSolution()
222 PetscCall(DMCreateGlobalVector(dm, &user->x)); in UpdateSolution()
230 if (param->ivisc == VISC_CONST && !param->stop_solve) { in UpdateSolution()
231 param->ivisc = VISC_CONST; in UpdateSolution()
233 PetscCall(SNESSolve(snes, 0, user->x)); in UpdateSolution()
237 PetscCall(VecCopy(user->x, user->Xguess)); in UpdateSolution()
238 if (param->stop_solve) goto done; in UpdateSolution()
242 if (param->ivisc >= VISC_DIFN && !param->stop_solve) { in UpdateSolution()
246 for (param->continuation = 0.0;; param->continuation += cont_incr) { in UpdateSolution()
247 …Call(PetscPrintf(PETSC_COMM_WORLD, " Continuation parameter = %g\n", (double)param->continuation)); in UpdateSolution()
249 /* solve the non-linear system */ in UpdateSolution()
250 PetscCall(VecCopy(user->Xguess, user->x)); in UpdateSolution()
251 PetscCall(SNESSolve(snes, 0, user->x)); in UpdateSolution()
256 if (param->stop_solve) goto done; in UpdateSolution()
260 cont_incr = -PetscAbsReal(cont_incr) / 2.0; in UpdateSolution()
265 PetscCall(VecCopy(user->x, user->Xguess)); in UpdateSolution()
266 if (param->continuation >= 1.0) goto done; in UpdateSolution()
271 if (param->continuation + cont_incr > 1.0) cont_incr = 1.0 - param->continuation; in UpdateSolution()
276 …if (param->stop_solve && !q) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "USER SIGNAL: stopping solve.… in UpdateSolution()
308 Parameter *param = user->param; in HorizVelocity()
309 GridInfo *grid = user->grid; in HorizVelocity()
310 PetscScalar st, ct, th, c = param->c, d = param->d; in HorizVelocity()
313 x = (i - grid->jlid) * grid->dx; in HorizVelocity()
314 z = (j - grid->jlid - 0.5) * grid->dz; in HorizVelocity()
319 return ct * (c * th * st + d * (st + th * ct)) + st * (c * (st - th * ct) + d * th * st); in HorizVelocity()
325 Parameter *param = user->param; in VertVelocity()
326 GridInfo *grid = user->grid; in VertVelocity()
327 PetscScalar st, ct, th, c = param->c, d = param->d; in VertVelocity()
330 x = (i - grid->jlid - 0.5) * grid->dx; in VertVelocity()
331 z = (j - grid->jlid) * grid->dz; in VertVelocity()
336 return st * (c * th * st + d * (st + th * ct)) - ct * (c * (st - th * ct) + d * th * st); in VertVelocity()
342 Parameter *param = user->param; in Pressure()
343 GridInfo *grid = user->grid; in Pressure()
344 PetscScalar x, z, r, st, ct, c = param->c, d = param->d; in Pressure()
346 x = (i - grid->jlid - 0.5) * grid->dx; in Pressure()
347 z = (j - grid->jlid - 0.5) * grid->dz; in Pressure()
351 return -2.0 * (c * ct - d * st) / r; in Pressure()
357 Parameter *param = user->param; in CalcSecInv()
358 GridInfo *grid = user->grid; in CalcSecInv()
359 PetscInt ilim = grid->ni - 1, jlim = grid->nj - 1; in CalcSecInv()
364 if (i == ilim) i--; in CalcSecInv()
365 if (j == jlim) j--; in CalcSecInv()
368 if (j <= grid->jlid) return EPS_ZERO; in CalcSecInv()
371 uW = x[j][i - 1].u; in CalcSecInv()
373 wS = x[j - 1][i].w; in CalcSecInv()
374 wE = WInterp(x, i, j - 1); in CalcSecInv()
376 uN = param->cb; in CalcSecInv()
377 wW = param->sb; in CalcSecInv()
379 uN = UInterp(x, i - 1, j); in CalcSecInv()
380 wW = WInterp(x, i - 1, j - 1); in CalcSecInv()
383 if (j == grid->jlid + 1) uS = 0.0; in CalcSecInv()
384 else uS = UInterp(x, i - 1, j - 1); in CalcSecInv()
387 if (j < grid->jlid) return EPS_ZERO; in CalcSecInv()
394 wN = param->sb; in CalcSecInv()
395 uW = param->cb; in CalcSecInv()
398 uW = UInterp(x, i - 1, j); in CalcSecInv()
401 if (j == grid->jlid) { in CalcSecInv()
404 uS = -uN; in CalcSecInv()
405 wS = -wN; in CalcSecInv()
408 wS = WInterp(x, i, j - 1); in CalcSecInv()
412 eps11 = (uE - uW) / grid->dx; in CalcSecInv()
413 eps22 = (wN - wS) / grid->dz; in CalcSecInv()
414 eps12 = 0.5 * ((uN - uS) / grid->dz + (wE - wW) / grid->dx); in CalcSecInv()
423 ViscParam difn = param->diffusion, disl = param->dislocation; in Viscosity()
424 PetscInt iVisc = param->ivisc; in Viscosity()
425 PetscScalar eps_scale = param->V / (param->L * 1000.0); in Viscosity()
429 P = rho_g * (z * param->L * 1000.0); /* Pa */ in Viscosity()
436 …scRealPart(difn.A * PetscExpScalar((difn.Estar + P * difn.Vstar) / R / (T + 273.0)) / param->eta0); in Viscosity()
439 strain_power = PetscPowScalar(eps * eps_scale, (1.0 - disl.n) / disl.n); in Viscosity()
441 …cExpScalar((disl.Estar + P * disl.Vstar) / disl.n / R / (T + 273.0)) * strain_power / param->eta0); in Viscosity()
444 strain_power = PetscPowScalar(eps * eps_scale, (1.0 - disl.n) / disl.n); in Viscosity()
446 v1 = difn.A * PetscExpScalar((difn.Estar + P * difn.Vstar) / R / (T + 273.0)) / param->eta0; in Viscosity()
447 …scExpScalar((disl.Estar + P * disl.Vstar) / disl.n / R / (T + 273.0)) * strain_power / param->eta0; in Viscosity()
452 /* max viscosity is param->eta0 */ in Viscosity()
454 /* min viscosity is param->visc_cutoff */ in Viscosity()
455 result = PetscMax(result, param->visc_cutoff); in Viscosity()
457 result = PetscPowReal(result, param->continuation); in Viscosity()
461 /* computes the residual of the x-component of eqn (1) above */
464 Parameter *param = user->param; in XMomentumResidual()
465 GridInfo *grid = user->grid; in XMomentumResidual()
466 PetscScalar dx = grid->dx, dz = grid->dz; in XMomentumResidual()
470 PetscInt jlim = grid->nj - 1; in XMomentumResidual()
472 z_scale = param->z_scale; in XMomentumResidual()
474 if (param->ivisc == VISC_DIFN || param->ivisc >= VISC_DISL) { /* viscosity is T-dependent */ in XMomentumResidual()
475 TS = param->potentialT * TInterp(x, i, j - 1) * PetscExpScalar((j - 1.0) * dz * z_scale); in XMomentumResidual()
477 else TN = param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * z_scale); in XMomentumResidual()
478 TW = param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * z_scale); in XMomentumResidual()
479 TE = param->potentialT * x[j][i + 1].T * PetscExpScalar((j - 0.5) * dz * z_scale); in XMomentumResidual()
480 if (param->ivisc >= VISC_DISL) { /* olivine dislocation creep */ in XMomentumResidual()
482 epsS = CalcSecInv(x, i, j - 1, CELL_CORNER, user); in XMomentumResidual()
488 etaS = Viscosity(TS, epsS, dz * (j - 0.5), param); in XMomentumResidual()
492 dPdx = (x[j][i + 1].p - x[j][i].p) / dx; in XMomentumResidual()
493 if (j == jlim) dudzN = etaN * (x[j][i].w - x[j][i + 1].w) / dx; in XMomentumResidual()
494 else dudzN = etaN * (x[j + 1][i].u - x[j][i].u) / dz; in XMomentumResidual()
495 dudzS = etaS * (x[j][i].u - x[j - 1][i].u) / dz; in XMomentumResidual()
496 dudxE = etaE * (x[j][i + 1].u - x[j][i].u) / dx; in XMomentumResidual()
497 dudxW = etaW * (x[j][i].u - x[j][i - 1].u) / dx; in XMomentumResidual()
499 residual = -dPdx /* X-MOMENTUM EQUATION*/ in XMomentumResidual()
500 + (dudxE - dudxW) / dx + (dudzN - dudzS) / dz; in XMomentumResidual()
502 if (param->ivisc != VISC_CONST) { in XMomentumResidual()
503 dwdxN = etaN * (x[j][i + 1].w - x[j][i].w) / dx; in XMomentumResidual()
504 dwdxS = etaS * (x[j - 1][i + 1].w - x[j - 1][i].w) / dx; in XMomentumResidual()
506 residual += (dudxE - dudxW) / dx + (dwdxN - dwdxS) / dz; in XMomentumResidual()
512 /* computes the residual of the z-component of eqn (1) above */
515 Parameter *param = user->param; in ZMomentumResidual()
516 GridInfo *grid = user->grid; in ZMomentumResidual()
517 PetscScalar dx = grid->dx, dz = grid->dz; in ZMomentumResidual()
521 PetscInt ilim = grid->ni - 1; in ZMomentumResidual()
524 z_scale = param->z_scale; in ZMomentumResidual()
527 if (param->ivisc == VISC_DIFN || param->ivisc >= VISC_DISL) { /* viscosity is T-dependent */ in ZMomentumResidual()
528 TN = param->potentialT * x[j + 1][i].T * PetscExpScalar((j + 0.5) * dz * z_scale); in ZMomentumResidual()
529 TS = param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * z_scale); in ZMomentumResidual()
530 TW = param->potentialT * TInterp(x, i - 1, j) * PetscExpScalar(j * dz * z_scale); in ZMomentumResidual()
532 else TE = param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * z_scale); in ZMomentumResidual()
533 if (param->ivisc >= VISC_DISL) { /* olivine dislocation creep */ in ZMomentumResidual()
537 epsW = CalcSecInv(x, i - 1, j, CELL_CORNER, user); in ZMomentumResidual()
545 dPdz = (x[j + 1][i].p - x[j][i].p) / dz; in ZMomentumResidual()
546 dwdzN = etaN * (x[j + 1][i].w - x[j][i].w) / dz; in ZMomentumResidual()
547 dwdzS = etaS * (x[j][i].w - x[j - 1][i].w) / dz; in ZMomentumResidual()
548 if (i == ilim) dwdxE = etaE * (x[j][i].u - x[j + 1][i].u) / dz; in ZMomentumResidual()
549 else dwdxE = etaE * (x[j][i + 1].w - x[j][i].w) / dx; in ZMomentumResidual()
550 dwdxW = 2.0 * etaW * (x[j][i].w - x[j][i - 1].w) / dx; in ZMomentumResidual()
552 /* Z-MOMENTUM */ in ZMomentumResidual()
553 residual = -dPdz /* constant viscosity terms */ in ZMomentumResidual()
554 + (dwdzN - dwdzS) / dz + (dwdxE - dwdxW) / dx; in ZMomentumResidual()
556 if (param->ivisc != VISC_CONST) { in ZMomentumResidual()
557 dudzE = etaE * (x[j + 1][i].u - x[j][i].u) / dz; in ZMomentumResidual()
558 dudzW = etaW * (x[j + 1][i - 1].u - x[j][i - 1].u) / dz; in ZMomentumResidual()
560 residual += (dwdzN - dwdzS) / dz + (dudzE - dudzW) / dx; in ZMomentumResidual()
569 GridInfo *grid = user->grid; in ContinuityResidual()
572 uW = x[j][i - 1].u; in ContinuityResidual()
574 dudx = (uE - uW) / grid->dx; in ContinuityResidual()
575 wS = x[j - 1][i].w; in ContinuityResidual()
577 dwdz = (wN - wS) / grid->dz; in ContinuityResidual()
585 Parameter *param = user->param; in EnergyResidual()
586 GridInfo *grid = user->grid; in EnergyResidual()
587 PetscScalar dx = grid->dx, dz = grid->dz; in EnergyResidual()
588 PetscInt ilim = grid->ni - 1, jlim = grid->nj - 1, jlid = grid->jlid; in EnergyResidual()
593 dTdzN = (x[j + 1][i].T - x[j][i].T) / dz; in EnergyResidual()
594 dTdzS = (x[j][i].T - x[j - 1][i].T) / dz; in EnergyResidual()
595 dTdxE = (x[j][i + 1].T - x[j][i].T) / dx; in EnergyResidual()
596 dTdxW = (x[j][i].T - x[j][i - 1].T) / dx; in EnergyResidual()
598 residual = ((dTdzN - dTdzS) / dz + /* diffusion term */ in EnergyResidual()
599 (dTdxE - dTdxW) / dx) * in EnergyResidual()
600 dx * dz / param->peclet; in EnergyResidual()
607 uW = uE = param->cb; in EnergyResidual()
608 wS = wN = param->sb; in EnergyResidual()
611 uW = x[j][i - 1].u; in EnergyResidual()
613 wS = x[j - 1][i].w; in EnergyResidual()
617 if (param->adv_scheme == ADVECT_FV || i == ilim - 1 || j == jlim - 1 || i == 1 || j == 1) { in EnergyResidual()
619 TS = (x[j][i].T + x[j - 1][i].T) / 2.0; in EnergyResidual()
622 TW = (x[j][i].T + x[j][i - 1].T) / 2.0; in EnergyResidual()
630 …-x[j][i + 2].T + 5.0 * (x[j][i + 1].T + x[j][i].T) - x[j][i - 1].T) / 8.0 - PetscAbsScalar(uE) * (… in EnergyResidual()
631 …-x[j][i + 1].T + 5.0 * (x[j][i].T + x[j][i - 1].T) - x[j][i - 2].T) / 8.0 - PetscAbsScalar(uW) * (… in EnergyResidual()
632 …-x[j + 2][i].T + 5.0 * (x[j + 1][i].T + x[j][i].T) - x[j - 1][i].T) / 8.0 - PetscAbsScalar(wN) * (… in EnergyResidual()
633 …-x[j + 1][i].T + 5.0 * (x[j][i].T + x[j - 1][i].T) - x[j - 2][i].T) / 8.0 - PetscAbsScalar(wS) * (… in EnergyResidual()
636 residual -= (fE - fW + fN - fS); in EnergyResidual()
641 /* computes the shear stress---used on the boundaries */
644 Parameter *param = user->param; in ShearStress()
645 GridInfo *grid = user->grid; in ShearStress()
646 PetscInt ilim = grid->ni - 1, jlim = grid->nj - 1; in ShearStress()
649 if (j <= grid->jlid || i < j || i == ilim || j == jlim) return EPS_ZERO; in ShearStress()
653 wE = WInterp(x, i, j - 1); in ShearStress()
655 wW = param->sb; in ShearStress()
656 uN = param->cb; in ShearStress()
658 wW = WInterp(x, i - 1, j - 1); in ShearStress()
659 uN = UInterp(x, i - 1, j); in ShearStress()
661 if (j == grid->jlid + 1) uS = 0.0; in ShearStress()
662 else uS = UInterp(x, i - 1, j - 1); in ShearStress()
672 return (uN - uS) / grid->dz + (wE - wW) / grid->dx; in ShearStress()
675 /* computes the normal stress---used on the boundaries */
678 Parameter *param = user->param; in XNormalStress()
679 GridInfo *grid = user->grid; in XNormalStress()
680 PetscScalar dx = grid->dx, dz = grid->dz; in XNormalStress()
681 PetscInt ilim = grid->ni - 1, jlim = grid->nj - 1, ivisc; in XNormalStress()
683 if (i < j || j <= grid->jlid) return EPS_ZERO; in XNormalStress()
685 ivisc = param->ivisc; in XNormalStress()
686 z_scale = param->z_scale; in XNormalStress()
690 TC = param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * z_scale); in XNormalStress()
694 uW = x[j][i - 1].u; in XNormalStress()
701 TC = param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * z_scale); in XNormalStress()
705 if (i == j) uW = param->sb; in XNormalStress()
706 else uW = UInterp(x, i - 1, j); in XNormalStress()
711 return 2.0 * etaC * (uE - uW) / dx - pC; in XNormalStress()
714 /* computes the normal stress---used on the boundaries */
717 Parameter *param = user->param; in ZNormalStress()
718 GridInfo *grid = user->grid; in ZNormalStress()
719 PetscScalar dz = grid->dz; in ZNormalStress()
720 PetscInt ilim = grid->ni - 1, jlim = grid->nj - 1, ivisc; in ZNormalStress()
723 if (i < j || j <= grid->jlid) return EPS_ZERO; in ZNormalStress()
725 ivisc = param->ivisc; in ZNormalStress()
726 z_scale = param->z_scale; in ZNormalStress()
730 TC = param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * z_scale); in ZNormalStress()
734 wS = x[j - 1][i].w; in ZNormalStress()
740 TC = param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * z_scale); in ZNormalStress()
743 if (i == j) wN = param->sb; in ZNormalStress()
745 wS = WInterp(x, i, j - 1); in ZNormalStress()
749 return 2.0 * etaC * (wN - wS) / dz - pC; in ZNormalStress()
753 INITIALIZATION, POST-PROCESSING AND OUTPUT FUNCTIONS
761 PetscReal alpha_g_on_cp_units_inverse_km = 4.0e-5 * 9.8; in SetParams()
765 param->slab_dip = 45.0; in SetParams()
766 param->width = 320.0; /* km */ in SetParams()
767 param->depth = 300.0; /* km */ in SetParams()
768 param->lid_depth = 35.0; /* km */ in SetParams()
769 param->fault_depth = 35.0; /* km */ in SetParams()
771 PetscCall(PetscOptionsGetReal(NULL, NULL, "-slab_dip", ¶m->slab_dip, NULL)); in SetParams()
772 PetscCall(PetscOptionsGetReal(NULL, NULL, "-width", ¶m->width, NULL)); in SetParams()
773 PetscCall(PetscOptionsGetReal(NULL, NULL, "-depth", ¶m->depth, NULL)); in SetParams()
774 PetscCall(PetscOptionsGetReal(NULL, NULL, "-lid_depth", ¶m->lid_depth, NULL)); in SetParams()
775 PetscCall(PetscOptionsGetReal(NULL, NULL, "-fault_depth", ¶m->fault_depth, NULL)); in SetParams()
777 param->slab_dip = param->slab_dip * PETSC_PI / 180.0; /* radians */ in SetParams()
780 PetscCall(PetscOptionsGetInt(NULL, NULL, "-jfault", &grid->jfault, NULL)); in SetParams()
781 grid->ni = 82; in SetParams()
782 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ni", &grid->ni, NULL)); in SetParams()
784 grid->dx = param->width / ((PetscReal)(grid->ni - 2)); /* km */ in SetParams()
785 grid->dz = grid->dx * PetscTanReal(param->slab_dip); /* km */ in SetParams()
786 grid->nj = (PetscInt)(param->depth / grid->dz + 3.0); /* gridpoints*/ in SetParams()
787 param->depth = grid->dz * (grid->nj - 2); /* km */ in SetParams()
788 grid->inose = 0; /* gridpoints*/ in SetParams()
789 PetscCall(PetscOptionsGetInt(NULL, NULL, "-inose", &grid->inose, NULL)); in SetParams()
790 grid->bx = DM_BOUNDARY_NONE; in SetParams()
791 grid->by = DM_BOUNDARY_NONE; in SetParams()
792 grid->stencil = DMDA_STENCIL_BOX; in SetParams()
793 grid->dof = 4; in SetParams()
794 grid->stencil_width = 2; in SetParams()
795 grid->mglevels = 1; in SetParams()
798 param->pv_analytic = PETSC_FALSE; in SetParams()
799 param->ibound = BC_NOSTRESS; in SetParams()
800 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ibound", ¶m->ibound, NULL)); in SetParams()
803 param->slab_velocity = 5.0; /* cm/yr */ in SetParams()
804 param->slab_age = 50.0; /* Ma */ in SetParams()
805 param->lid_age = 50.0; /* Ma */ in SetParams()
806 param->kappa = 0.7272e-6; /* m^2/sec */ in SetParams()
807 param->potentialT = 1300.0; /* degrees C */ in SetParams()
809 PetscCall(PetscOptionsGetReal(NULL, NULL, "-slab_velocity", ¶m->slab_velocity, NULL)); in SetParams()
810 PetscCall(PetscOptionsGetReal(NULL, NULL, "-slab_age", ¶m->slab_age, NULL)); in SetParams()
811 PetscCall(PetscOptionsGetReal(NULL, NULL, "-lid_age", ¶m->lid_age, NULL)); in SetParams()
812 PetscCall(PetscOptionsGetReal(NULL, NULL, "-kappa", ¶m->kappa, NULL)); in SetParams()
813 PetscCall(PetscOptionsGetReal(NULL, NULL, "-potentialT", ¶m->potentialT, NULL)); in SetParams()
816 param->ivisc = 3; /* 0=isovisc, 1=difn creep, 2=disl creep, 3=full */ in SetParams()
817 param->eta0 = 1e24; /* Pa-s */ in SetParams()
818 param->visc_cutoff = 0.0; /* factor of eta_0 */ in SetParams()
819 param->continuation = 1.0; in SetParams()
822 param->diffusion.A = 1.8e7; /* Pa-s */ in SetParams()
823 param->diffusion.n = 1.0; /* dim'less */ in SetParams()
824 param->diffusion.Estar = 375e3; /* J/mol */ in SetParams()
825 param->diffusion.Vstar = 5e-6; /* m^3/mol */ in SetParams()
827 /* constants for param->dislocationocation creep */ in SetParams()
828 param->dislocation.A = 2.8969e4; /* Pa-s */ in SetParams()
829 param->dislocation.n = 3.5; /* dim'less */ in SetParams()
830 param->dislocation.Estar = 530e3; /* J/mol */ in SetParams()
831 param->dislocation.Vstar = 14e-6; /* m^3/mol */ in SetParams()
833 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ivisc", ¶m->ivisc, NULL)); in SetParams()
834 PetscCall(PetscOptionsGetReal(NULL, NULL, "-visc_cutoff", ¶m->visc_cutoff, NULL)); in SetParams()
836 param->output_ivisc = param->ivisc; in SetParams()
838 PetscCall(PetscOptionsGetInt(NULL, NULL, "-output_ivisc", ¶m->output_ivisc, NULL)); in SetParams()
839 PetscCall(PetscOptionsGetReal(NULL, NULL, "-vstar", ¶m->dislocation.Vstar, NULL)); in SetParams()
842 param->quiet = PETSC_FALSE; in SetParams()
843 param->param_test = PETSC_FALSE; in SetParams()
845 PetscCall(PetscOptionsHasName(NULL, NULL, "-quiet", ¶m->quiet)); in SetParams()
846 PetscCall(PetscOptionsHasName(NULL, NULL, "-test", ¶m->param_test)); in SetParams()
847 …PetscCall(PetscOptionsGetString(NULL, NULL, "-file", param->filename, sizeof(param->filename), &pa… in SetParams()
850 param->adv_scheme = ADVECT_FROMM; /* advection scheme: 0=finite vol, 1=Fromm */ in SetParams()
852 PetscCall(PetscOptionsGetInt(NULL, NULL, "-adv_scheme", ¶m->adv_scheme, NULL)); in SetParams()
855 param->stop_solve = PETSC_FALSE; in SetParams()
856 param->interrupted = PETSC_FALSE; in SetParams()
857 param->kspmon = PETSC_FALSE; in SetParams()
858 param->toggle_kspmon = PETSC_FALSE; in SetParams()
861 param->sb = PetscSinReal(param->slab_dip); in SetParams()
862 param->cb = PetscCosReal(param->slab_dip); in SetParams()
863 …param->c = param->slab_dip * param->sb / (param->slab_dip * param->slab_dip - param->sb * param->… in SetParams()
864 …param->d = (param->slab_dip * param->cb - param->sb) / (param->slab_dip * param->slab_dip - param… in SetParams()
866 /* length, velocity and time scale for non-dimensionalization */ in SetParams()
867 param->L = PetscMin(param->width, param->depth); /* km */ in SetParams()
868 param->V = param->slab_velocity / 100.0 / SEC_PER_YR; /* m/sec */ in SetParams()
871 param->scaled_width = param->width / param->L; /* dim'less */ in SetParams()
872 param->scaled_depth = param->depth / param->L; /* dim'less */ in SetParams()
873 param->lid_depth = param->lid_depth / param->L; /* dim'less */ in SetParams()
874 param->fault_depth = param->fault_depth / param->L; /* dim'less */ in SetParams()
875 grid->dx = grid->dx / param->L; /* dim'less */ in SetParams()
876 grid->dz = grid->dz / param->L; /* dim'less */ in SetParams()
877 grid->jlid = (PetscInt)(param->lid_depth / grid->dz); /* gridcells */ in SetParams()
878 grid->jfault = (PetscInt)(param->fault_depth / grid->dz); /* gridcells */ in SetParams()
879 param->lid_depth = grid->jlid * grid->dz; /* dim'less */ in SetParams()
880 param->fault_depth = grid->jfault * grid->dz; /* dim'less */ in SetParams()
881 grid->corner = grid->jlid + 1; /* gridcells */ in SetParams()
882 param->peclet = param->V /* m/sec */ in SetParams()
883 * param->L * 1000.0 /* m */ in SetParams()
884 / param->kappa; /* m^2/sec */ in SetParams()
885 param->z_scale = param->L * alpha_g_on_cp_units_inverse_km; in SetParams()
886 param->skt = PetscSqrtReal(param->kappa * param->slab_age * SEC_PER_YR); in SetParams()
887 PetscCall(PetscOptionsGetReal(NULL, NULL, "-peclet", ¶m->peclet, NULL)); in SetParams()
899 if (!param->quiet) { in ReportParams()
900 …PetscCall(PetscPrintf(PETSC_COMM_WORLD, "---------------------BEGIN ex30 PARAM REPORT-------------… in ReportParams()
902 …MM_WORLD, " Width = %g km, Depth = %g km\n", (double)param->width, (double)param->depth)); in ReportParams()
903 …degrees, Slab velocity = %g cm/yr\n", (double)(param->slab_dip * 180.0 / PETSC_PI), (double)param… in ReportParams()
904 …m, Fault depth = %5.2f km\n", (double)(param->lid_depth * param->L), (double)(param->fault_depth… in ReportParams()
907 …FMT " [dx,dz] = %g, %g km\n", grid->ni, grid->nj, (double)(grid->dx * param->L), (double)(gr… in ReportParams()
908 …jlid = %3" PetscInt_FMT " jfault = %3" PetscInt_FMT " \n", grid->jlid, grid->jfault)); in ReportParams()
909 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Pe = %g\n", (double)param->peclet)); in ReportParams()
912 if (param->ivisc == VISC_CONST) { in ReportParams()
914 …if (param->pv_analytic) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Pressur… in ReportParams()
915 } else if (param->ivisc == VISC_DIFN) { in ReportParams()
916 …PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Diffusion Creep (T-Dependent Newtonian) … in ReportParams()
917 … Viscosity range: %g--%g Pa-sec \n", (double)param->eta0, (double)(param->visc_… in ReportParams()
918 } else if (param->ivisc == VISC_DISL) { in ReportParams()
919 …PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Dislocation Creep (T-Dependent Non-Newto… in ReportParams()
920 … Viscosity range: %g--%g Pa-sec \n", (double)param->eta0, (double)(param->visc_… in ReportParams()
921 } else if (param->ivisc == VISC_FULL) { in ReportParams()
923 … Viscosity range: %g--%g Pa-sec \n", (double)param->eta0, (double)(param->visc_… in ReportParams()
930 if (param->ibound == BC_ANALYTIC) { in ReportParams()
932 } else if (param->ibound == BC_NOSTRESS) { in ReportParams()
933 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Stress-Free (normal & shear stress)\n")); in ReportParams()
934 } else if (param->ibound == BC_EXPERMNT) { in ReportParams()
935 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Experimental boundary condition \n")); in ReportParams()
941 if (param->output_to_file) { in ReportParams()
943 …PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Output Destination: Mat file \"%s\"\n", param->file… in ReportParams()
945 …Printf(PETSC_COMM_WORLD, "Output Destination: PETSc binary file \"%s\"\n", param->filename)); in ReportParams()
948 …->output_ivisc != param->ivisc) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " … in ReportParams()
950 …PetscCall(PetscPrintf(PETSC_COMM_WORLD, "---------------------END ex30 PARAM REPORT---------------… in ReportParams()
952 if (param->param_test) PetscCall(PetscEnd()); in ReportParams()
956 /* ------------------------------------------------------------------- */
960 /* ------------------------------------------------------------------- */ in Initialize()
972 Xguess = user->Xguess; in Initialize()
973 param = user->param; in Initialize()
974 grid = user->grid; in Initialize()
981 if (i < j) x[j][i].u = param->cb; in Initialize()
982 else if (j <= grid->jlid) x[j][i].u = 0.0; in Initialize()
985 if (i <= j) x[j][i].w = param->sb; in Initialize()
986 else if (j <= grid->jlid) x[j][i].w = 0.0; in Initialize()
989 if (i < j || j <= grid->jlid) x[j][i].p = 0.0; in Initialize()
992 x[j][i].T = PetscMin(grid->dz * (j - 0.5), 1.0); in Initialize()
1017 param = user->param; in DoOutput()
1018 grid = user->grid; in DoOutput()
1019 ivt = param->ivisc; in DoOutput()
1021 param->ivisc = param->output_ivisc; in DoOutput()
1025 PetscCall(ViscosityField(da, user->x, user->Xguess)); in DoOutput()
1031 if (param->output_to_file) { /* send output to binary file */ in DoOutput()
1036 PetscCall(VecSetValue(pars, 0, (PetscScalar)grid->ni, INSERT_VALUES)); in DoOutput()
1037 PetscCall(VecSetValue(pars, 1, (PetscScalar)grid->nj, INSERT_VALUES)); in DoOutput()
1038 PetscCall(VecSetValue(pars, 2, (PetscScalar)grid->dx, INSERT_VALUES)); in DoOutput()
1039 PetscCall(VecSetValue(pars, 3, (PetscScalar)grid->dz, INSERT_VALUES)); in DoOutput()
1040 PetscCall(VecSetValue(pars, 4, (PetscScalar)param->L, INSERT_VALUES)); in DoOutput()
1041 PetscCall(VecSetValue(pars, 5, (PetscScalar)param->V, INSERT_VALUES)); in DoOutput()
1043 PetscCall(VecSetValue(pars, 7, (PetscScalar)param->slab_dip, INSERT_VALUES)); in DoOutput()
1044 PetscCall(VecSetValue(pars, 8, (PetscScalar)grid->jlid, INSERT_VALUES)); in DoOutput()
1045 PetscCall(VecSetValue(pars, 9, (PetscScalar)param->lid_depth, INSERT_VALUES)); in DoOutput()
1046 PetscCall(VecSetValue(pars, 10, (PetscScalar)grid->jfault, INSERT_VALUES)); in DoOutput()
1047 PetscCall(VecSetValue(pars, 11, (PetscScalar)param->fault_depth, INSERT_VALUES)); in DoOutput()
1048 PetscCall(VecSetValue(pars, 12, (PetscScalar)param->potentialT, INSERT_VALUES)); in DoOutput()
1049 PetscCall(VecSetValue(pars, 13, (PetscScalar)param->ivisc, INSERT_VALUES)); in DoOutput()
1050 PetscCall(VecSetValue(pars, 14, (PetscScalar)param->visc_cutoff, INSERT_VALUES)); in DoOutput()
1051 PetscCall(VecSetValue(pars, 15, (PetscScalar)param->ibound, INSERT_VALUES)); in DoOutput()
1062 PetscCall(PetscViewerMatlabOpen(PETSC_COMM_WORLD, param->filename, FILE_MODE_WRITE, &viewer)); in DoOutput()
1064 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, param->filename, FILE_MODE_WRITE, &viewer)); in DoOutput()
1070 PetscCall(PetscObjectSetName((PetscObject)user->x, "out")); in DoOutput()
1071 PetscCall(VecView(user->x, viewer)); in DoOutput()
1072 PetscCall(PetscObjectSetName((PetscObject)user->Xguess, "aux")); in DoOutput()
1073 PetscCall(VecView(user->Xguess, viewer)); in DoOutput()
1075 PetscCall(PetscObjectSetName((PetscObject)user->Xguess, "str")); in DoOutput()
1076 PetscCall(VecView(user->Xguess, viewer)); in DoOutput()
1085 param->ivisc = ivt; in DoOutput()
1089 /* ------------------------------------------------------------------- */
1092 /* ------------------------------------------------------------------- */ in ViscosityField()
1104 param = user->param; in ViscosityField()
1105 grid = user->grid; in ViscosityField()
1106 ivt = param->ivisc; in ViscosityField()
1107 param->ivisc = param->output_ivisc; in ViscosityField()
1116 /* dx = grid->dx; */ dz = grid->dz; in ViscosityField()
1118 ilim = grid->ni - 1; in ViscosityField()
1119 jlim = grid->nj - 1; in ViscosityField()
1125 …T = PetscRealPart(param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * param->z_scale)); in ViscosityField()
1127 …TC = PetscRealPart(param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * param->z_scale)); in ViscosityField()
1136 v[j][i].p = Viscosity(T, eps, dz * (j - 0.5), param); in ViscosityField()
1144 param->ivisc = ivt; in ViscosityField()
1148 /* ------------------------------------------------------------------- */
1149 /* post-processing: compute stress everywhere */
1151 /* ------------------------------------------------------------------- */ in StressField()
1163 PetscCall(DMDAVecGetArray(da, user->Xguess, (void **)&x)); in StressField()
1166 PetscCall(DMGlobalToLocalBegin(da, user->x, INSERT_VALUES, locVec)); in StressField()
1167 PetscCall(DMGlobalToLocalEnd(da, user->x, INSERT_VALUES, locVec)); in StressField()
1181 PetscCall(DMDAVecRestoreArray(da, user->Xguess, (void **)&x)); in StressField()
1194 Parameter *param = user->param; in SlabVel()
1195 GridInfo *grid = user->grid; in SlabVel()
1198 if (i < j - 1) return param->cb; in SlabVel()
1199 else if (j <= grid->jfault) return 0.0; in SlabVel()
1200 else return param->cb; in SlabVel()
1203 if (i < j) return param->sb; in SlabVel()
1204 else if (j <= grid->jfault) return 0.0; in SlabVel()
1205 else return param->sb; in SlabVel()
1209 /* solution to diffusive half-space cooling model for BC */
1212 Parameter *param = user->param; in PlateModel()
1214 if (plate == PLATE_LID) z = (j - 0.5) * user->grid->dz; in PlateModel()
1215 else z = (j - 0.5) * user->grid->dz * param->cb; /* PLATE_SLAB */ in PlateModel()
1217 return (PetscReal)(erf((double)PetscRealPart(z * param->L / 2.0 / param->skt))); in PlateModel()
1228 /* ------------------------------------------------------------------- */
1230 /* ------------------------------------------------------------------- */ in SNESConverged_Interactive()
1233 Parameter *param = user->param; in SNESConverged_Interactive()
1237 if (param->interrupted) { in SNESConverged_Interactive()
1238 param->interrupted = PETSC_FALSE; in SNESConverged_Interactive()
1242 } else if (param->toggle_kspmon) { in SNESConverged_Interactive()
1243 param->toggle_kspmon = PETSC_FALSE; in SNESConverged_Interactive()
1247 if (param->kspmon) { in SNESConverged_Interactive()
1250 param->kspmon = PETSC_FALSE; in SNESConverged_Interactive()
1257 param->kspmon = PETSC_TRUE; in SNESConverged_Interactive()
1265 /* ------------------------------------------------------------------- */
1268 /* ------------------------------------------------------------------- */ in InteractiveHandler()
1271 Parameter *param = user->param; in InteractiveHandler()
1274 param->toggle_kspmon = PETSC_TRUE; in InteractiveHandler()
1277 param->interrupted = PETSC_TRUE; in InteractiveHandler()
1281 param->stop_solve = PETSC_TRUE; in InteractiveHandler()
1287 /* main call-back function that computes the processor-local piece of the residual */
1291 Parameter *param = user->param; in FormFunctionLocal()
1292 GridInfo *grid = user->grid; in FormFunctionLocal()
1299 mx = info->mx; in FormFunctionLocal()
1300 mz = info->my; in FormFunctionLocal()
1301 ilim = mx - 1; in FormFunctionLocal()
1302 jlim = mz - 1; in FormFunctionLocal()
1303 is = info->xs; in FormFunctionLocal()
1304 ie = info->xs + info->xm; in FormFunctionLocal()
1305 js = info->ys; in FormFunctionLocal()
1306 je = info->ys + info->ym; in FormFunctionLocal()
1309 /* ivisc = param->ivisc; */ ibound = param->ibound; in FormFunctionLocal()
1313 /************* X-MOMENTUM/VELOCITY *************/ in FormFunctionLocal()
1314 if (i < j) f[j][i].u = x[j][i].u - SlabVel('U', i, j, user); in FormFunctionLocal()
1315 … else if (j <= grid->jlid || (j < grid->corner + grid->inose && i < grid->corner + grid->inose)) { in FormFunctionLocal()
1317 f[j][i].u = x[j][i].u - 0.0; in FormFunctionLocal()
1321 f[j][i].u = x[j][i].u - HorizVelocity(i, j, user); in FormFunctionLocal()
1323 f[j][i].u = XNormalStress(x, i, j, CELL_CENTER, user) - EPS_ZERO; in FormFunctionLocal()
1329 f[j][i].u = x[j][i].u - HorizVelocity(i, j, user); in FormFunctionLocal()
1333 /* experimental boundary condition */ in FormFunctionLocal()
1341 /************* Z-MOMENTUM/VELOCITY *************/ in FormFunctionLocal()
1343 f[j][i].w = x[j][i].w - SlabVel('W', i, j, user); in FormFunctionLocal()
1345 …} else if (j <= grid->jlid || (j < grid->corner + grid->inose && i < grid->corner + grid->inose)) { in FormFunctionLocal()
1347 f[j][i].w = x[j][i].w - 0.0; in FormFunctionLocal()
1352 f[j][i].w = x[j][i].w - VertVelocity(i, j, user); in FormFunctionLocal()
1354 f[j][i].w = ZNormalStress(x, i, j, CELL_CENTER, user) - EPS_ZERO; in FormFunctionLocal()
1360 f[j][i].w = x[j][i].w - VertVelocity(i, j, user); in FormFunctionLocal()
1364 /* experimental boundary condition */ in FormFunctionLocal()
1373 …if (i < j || j <= grid->jlid || (j < grid->corner + grid->inose && i < grid->corner + grid->inose)… in FormFunctionLocal()
1379 f[j][i].p = x[j][i].p - Pressure(i, j, user); in FormFunctionLocal()
1393 f[j][i].T = x[j][i].T - PlateModel(j, PLATE_SLAB, user); in FormFunctionLocal()
1397 …mag_u = 1.0 - PetscPowRealInt(1.0 - PetscMax(PetscMin(PetscRealPart(x[j][i - 1].u) / param->cb… in FormFunctionLocal()
1398 …f[j][i].T = x[j][i].T - mag_u * x[j - 1][i - 1].T - (1.0 - mag_u) * PlateModel(j, PLATE_LID, user); in FormFunctionLocal()
1402 …mag_w = 1.0 - PetscPowRealInt(1.0 - PetscMax(PetscMin(PetscRealPart(x[j - 1][i].w) / param->sb… in FormFunctionLocal()
1403 f[j][i].T = x[j][i].T - mag_w * x[j - 1][i - 1].T - (1.0 - mag_w); in FormFunctionLocal()
1414 /*TEST
1419 test:
1420 args: -ni 18 -fp_trap 0
1421 filter: grep -v Destination
1424 TEST*/