Lines Matching refs:thi
299 static void THIInitialize_HOM_A(THI thi, PetscReal x, PetscReal y, PrmNode *p) in THIInitialize_HOM_A() argument
301 Units units = thi->units; in THIInitialize_HOM_A()
302 PetscReal s = -x * PetscSinReal(thi->alpha); in THIInitialize_HOM_A()
304 …500 * units->meter * PetscSinReal(x * 2 * PETSC_PI / thi->Lx) * PetscSinReal(y * 2 * PETSC_PI / th… in THIInitialize_HOM_A()
309 static void THIInitialize_HOM_C(THI thi, PetscReal x, PetscReal y, PrmNode *p) in THIInitialize_HOM_C() argument
311 Units units = thi->units; in THIInitialize_HOM_C()
312 PetscReal s = -x * PetscSinReal(thi->alpha); in THIInitialize_HOM_C()
317 …>beta2 = 1000 * (1 + PetscSinReal(x * 2 * PETSC_PI / thi->Lx) * PetscSinReal(y * 2 * PETSC_PI / th… in THIInitialize_HOM_C()
323 static void THIInitialize_HOM_X(THI thi, PetscReal xx, PetscReal yy, PrmNode *p) in THIInitialize_HOM_X() argument
325 Units units = thi->units; in THIInitialize_HOM_X()
326 …PetscReal x = xx * 2 * PETSC_PI / thi->Lx - PETSC_PI, y = yy * 2 * PETSC_PI / thi->Ly - PETSC_PI; … in THIInitialize_HOM_X()
327 PetscReal r = PetscSqrtReal(x * x + y * y), s = -x * PetscSinReal(thi->alpha); in THIInitialize_HOM_X()
334 static void THIInitialize_HOM_Y(THI thi, PetscReal xx, PetscReal yy, PrmNode *p) in THIInitialize_HOM_Y() argument
336 Units units = thi->units; in THIInitialize_HOM_Y()
337 …PetscReal x = xx * 2 * PETSC_PI / thi->Lx - PETSC_PI, y = yy * 2 * PETSC_PI / thi->Ly - PETSC_PI; … in THIInitialize_HOM_Y()
338 PetscReal r = PetscSqrtReal(x * x + y * y), s = -x * PetscSinReal(thi->alpha); in THIInitialize_HOM_Y()
347 static void THIInitialize_HOM_Z(THI thi, PetscReal xx, PetscReal yy, PrmNode *p) in THIInitialize_HOM_Z() argument
349 Units units = thi->units; in THIInitialize_HOM_Z()
350 …PetscReal x = xx * 2 * PETSC_PI / thi->Lx - PETSC_PI, y = yy * 2 * PETSC_PI / thi->Ly - PETSC_PI; … in THIInitialize_HOM_Z()
351 PetscReal r = PetscSqrtReal(x * x + y * y), s = -x * PetscSinReal(thi->alpha); in THIInitialize_HOM_Z()
358 static void THIFriction(THI thi, PetscReal rbeta2, PetscReal gam, PetscReal *beta2, PetscReal *dbet… in THIFriction() argument
360 if (thi->friction.irefgam == 0) { in THIFriction()
361 Units units = thi->units; in THIFriction()
362 … thi->friction.irefgam = 1. / (0.5 * PetscSqr(thi->friction.refvel * units->meter / units->year)); in THIFriction()
363 …thi->friction.eps2 = 0.5 * PetscSqr(thi->friction.epsvel * units->meter / units->year) * thi->f… in THIFriction()
365 if (thi->friction.exponent == 0) { in THIFriction()
369 …*beta2 = rbeta2 * PetscPowReal(thi->friction.eps2 + gam * thi->friction.irefgam, thi->friction.ex… in THIFriction()
370 …*dbeta2 = thi->friction.exponent * *beta2 / (thi->friction.eps2 + gam * thi->friction.irefgam) * t… in THIFriction()
374 static void THIViscosity(THI thi, PetscReal gam, PetscReal *eta, PetscReal *deta) in THIViscosity() argument
377 if (thi->viscosity.Bd2 == 0) { in THIViscosity()
378 Units units = thi->units; in THIViscosity()
383 thi->viscosity.Bd2 = B / 2; in THIViscosity()
384 thi->viscosity.exponent = (p - 2) / 2; in THIViscosity()
385 thi->viscosity.eps = 0.5 * PetscSqr(1e-5 / units->year); in THIViscosity()
387 Bd2 = thi->viscosity.Bd2; in THIViscosity()
388 exponent = thi->viscosity.exponent; in THIViscosity()
389 eps = thi->viscosity.eps; in THIViscosity()
416 static PetscErrorCode THIDestroy(THI *thi) in THIDestroy() argument
419 if (!*thi) PetscFunctionReturn(PETSC_SUCCESS); in THIDestroy()
420 if (--((PetscObject)*thi)->refct > 0) { in THIDestroy()
421 *thi = 0; in THIDestroy()
424 PetscCall(PetscFree((*thi)->units)); in THIDestroy()
425 PetscCall(PetscFree((*thi)->mattype)); in THIDestroy()
426 PetscCall(PetscHeaderDestroy(thi)); in THIDestroy()
433 THI thi; in THICreate() local
442 …PetscCall(PetscHeaderCreate(thi, THI_CLASSID, "THI", "Toy Hydrostatic Ice", "", comm, THIDestroy, … in THICreate()
444 PetscCall(PetscNew(&thi->units)); in THICreate()
445 units = thi->units; in THICreate()
460 thi->Lx = 10.e3; in THICreate()
461 thi->Ly = 10.e3; in THICreate()
462 thi->Lz = 1000; in THICreate()
463 thi->dirichlet_scale = 1; in THICreate()
464 thi->verbose = PETSC_FALSE; in THICreate()
473 L = thi->Lx; in THICreate()
475 if (flg) thi->Lx = thi->Ly = L; in THICreate()
476 PetscCall(PetscOptionsReal("-thi_Lx", "X Domain size (m)", "", thi->Lx, &thi->Lx, NULL)); in THICreate()
477 PetscCall(PetscOptionsReal("-thi_Ly", "Y Domain size (m)", "", thi->Ly, &thi->Ly, NULL)); in THICreate()
478 PetscCall(PetscOptionsReal("-thi_Lz", "Z Domain size (m)", "", thi->Lz, &thi->Lz, NULL)); in THICreate()
482 thi->initialize = THIInitialize_HOM_A; in THICreate()
483 thi->no_slip = PETSC_TRUE; in THICreate()
484 thi->alpha = 0.5; in THICreate()
487 thi->initialize = THIInitialize_HOM_C; in THICreate()
488 thi->no_slip = PETSC_FALSE; in THICreate()
489 thi->alpha = 0.1; in THICreate()
492 thi->initialize = THIInitialize_HOM_X; in THICreate()
493 thi->no_slip = PETSC_FALSE; in THICreate()
494 thi->alpha = 0.3; in THICreate()
497 thi->initialize = THIInitialize_HOM_Y; in THICreate()
498 thi->no_slip = PETSC_FALSE; in THICreate()
499 thi->alpha = 0.5; in THICreate()
502 thi->initialize = THIInitialize_HOM_Z; in THICreate()
503 thi->no_slip = PETSC_FALSE; in THICreate()
504 thi->alpha = 0.5; in THICreate()
520 …PetscCall(PetscOptionsReal("-thi_alpha", "Bed angle (degrees)", "", thi->alpha, &thi->alpha, NULL)… in THICreate()
522 thi->friction.refvel = 100.; in THICreate()
523 thi->friction.epsvel = 1.; in THICreate()
525 …("-thi_friction_refvel", "Reference velocity for sliding", "", thi->friction.refvel, &thi->frictio… in THICreate()
526 …i_friction_epsvel", "Regularization velocity for sliding", "", thi->friction.epsvel, &thi->frictio… in THICreate()
529 thi->friction.exponent = (m - 1) / 2; in THICreate()
531 …le", "Scale Dirichlet boundary conditions by this factor", "", thi->dirichlet_scale, &thi->dirichl… in THICreate()
532 …boundary conditions by this factor in SSA (2D) assembly", "", thi->ssa_friction_scale, &thi->ssa_f… in THICreate()
533 …"-thi_coarse2d", "Use a 2D coarse space corresponding to SSA", "", thi->coarse2d, &thi->coarse2d, … in THICreate()
534 …diagonal system (column coupling only) on the finest level", "", thi->tridiagonal, &thi->tridiagon… in THICreate()
536 PetscCall(PetscStrallocpy(mtype, (char **)&thi->mattype)); in THICreate()
537 …", "Enable verbose output (like matrix sizes and statistics)", "", thi->verbose, &thi->verbose, NU… in THICreate()
542 thi->Lx *= units->meter; in THICreate()
543 thi->Ly *= units->meter; in THICreate()
544 thi->Lz *= units->meter; in THICreate()
545 thi->alpha *= PETSC_PI / 180; in THICreate()
547 PRangeClear(&thi->eta); in THICreate()
548 PRangeClear(&thi->beta2); in THICreate()
552 driving = rho * grav * PetscSinReal(thi->alpha) * 1000 * units->meter; in THICreate()
553 THIViscosity(thi, 0.5 * gradu * gradu, &eta, &deta); in THICreate()
554 thi->rhog = rho * grav; in THICreate()
555 if (thi->verbose) { in THICreate()
556 …PetscCall(PetscPrintf(PetscObjectComm((PetscObject)thi), "Units: meter %8.2g second %8.2g kg %8.… in THICreate()
557 …mm((PetscObject)thi), "Domain (%6.2g,%6.2g,%6.2g), pressure %8.2g, driving stress %8.2g\n", (doubl… in THICreate()
558 …PetscCall(PetscPrintf(PetscObjectComm((PetscObject)thi), "Large velocity 1km/a %8.2g, velocity gra… in THICreate()
559 THIViscosity(thi, 0.5 * PetscSqr(1e-3 * gradu), &eta, &deta); in THICreate()
560 …PetscCall(PetscPrintf(PetscObjectComm((PetscObject)thi), "Small velocity 1m/a %8.2g, velocity gra… in THICreate()
564 *inthi = thi; in THICreate()
568 static PetscErrorCode THIInitializePrm(THI thi, DM da2prm, Vec prm) in THIInitializePrm() argument
579 PetscReal xx = thi->Lx * i / mx, yy = thi->Ly * j / my; in THIInitializePrm()
580 thi->initialize(thi, xx, yy, &p[i][j]); in THIInitializePrm()
587 static PetscErrorCode THISetUpDM(THI thi, DM dm) in THISetUpDM() argument
600 …PetscCall(DMDACreate2d(PetscObjectComm((PetscObject)thi), DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIOD… in THISetUpDM()
604 …PetscReal Lx = thi->Lx / thi->units->meter, Ly = thi->Ly / thi->units->meter, Lz = thi->Lz / thi->… in THISetUpDM()
606 …PetscCall(PetscPrintf(PetscObjectComm((PetscObject)thi), "Level %" PetscInt_FMT " domain size (m) … in THISetUpDM()
608 …PetscCall(PetscPrintf(PetscObjectComm((PetscObject)thi), "Level %" PetscInt_FMT " domain size (m) … in THISetUpDM()
611 PetscCall(THIInitializePrm(thi, da2prm, X)); in THISetUpDM()
612 if (thi->tridiagonal) { /* Reset coarse Jacobian evaluation */ in THISetUpDM()
613 PetscCall(DMDASNESSetJacobianLocal(dm, (DMDASNESJacobianFn *)THIJacobianLocal_3D_Full, thi)); in THISetUpDM()
615 …if (thi->coarse2d) PetscCall(DMDASNESSetJacobianLocal(dm, (DMDASNESJacobianFn *)THIJacobianLocal_2… in THISetUpDM()
625 THI thi = (THI)ctx; in DMCoarsenHook_THI() local
629 PetscCall(THISetUpDM(thi, dmc)); in DMCoarsenHook_THI()
633 PetscCall(DMCoarsenHookAdd(dmc, DMCoarsenHook_THI, NULL, thi)); in DMCoarsenHook_THI()
639 THI thi = (THI)ctx; in DMRefineHook_THI() local
642 PetscCall(THISetUpDM(thi, dmf)); in DMRefineHook_THI()
643 PetscCall(DMSetMatType(dmf, thi->mattype)); in DMRefineHook_THI()
644 PetscCall(DMRefineHookAdd(dmf, DMRefineHook_THI, NULL, thi)); in DMRefineHook_THI()
646 PetscCall(DMCoarsenHookAdd(dmf, DMCoarsenHook_THI, NULL, thi)); in DMRefineHook_THI()
680 THI thi; in THIInitial() local
689 PetscCall(DMGetApplicationContext(da, &thi)); in THIInitial()
694 hx = thi->Lx / mx; in THIInitial()
695 hy = thi->Ly / my; in THIInitial()
699 …= zm - 1, drivingx = thi->rhog * (prm[i + 1][j].b + prm[i + 1][j].h - prm[i - 1][j].b - prm[i - 1]… in THIInitial()
710 static void PointwiseNonlinearity(THI thi, const Node n[PETSC_RESTRICT], const PetscReal phi[PETSC_… in PointwiseNonlinearity() argument
728 THIViscosity(thi, PetscRealPart(gam), eta, deta); in PointwiseNonlinearity()
731 static void PointwiseNonlinearity2D(THI thi, Node n[], PetscReal phi[], PetscReal dphi[4][2], Petsc… in PointwiseNonlinearity2D() argument
749 THIViscosity(thi, PetscRealPart(gam), eta, deta); in PointwiseNonlinearity2D()
752 static PetscErrorCode THIFunctionLocal(DMDALocalInfo *info, Node ***x, Node ***f, THI thi) in THIFunctionLocal() argument
764 hx = thi->Lx / info->mz; in THIFunctionLocal()
765 hy = thi->Ly / info->my; in THIFunctionLocal()
785 if (thi->no_slip && k == 0) { in THIFunctionLocal()
795 PointwiseNonlinearity(thi, n, phi, dphi, &u, &v, du, dv, &eta, &deta); in THIFunctionLocal()
796 jw /= thi->rhog; /* scales residuals to be O(1) */ in THIFunctionLocal()
800 const PetscReal ds[2] = {-PetscSinReal(thi->alpha), 0}; in THIFunctionLocal()
802 …[1]) + dp[1] * jw * eta * (du[1] + dv[0]) + dp[2] * jw * eta * du[2] + pp * jw * thi->rhog * ds[0]; in THIFunctionLocal()
803 …[1]) + dp[0] * jw * eta * (du[1] + dv[0]) + dp[2] * jw * eta * dv[2] + pp * jw * thi->rhog * ds[1]; in THIFunctionLocal()
807 if (thi->no_slip) { in THIFunctionLocal()
822 …etscScalar diagu = 2 * etabase / thi->rhog * (hx * hy / hz + hx * hz / hy + 4 * hy * hz / hx), dia… in THIFunctionLocal()
823 fn[0]->u = thi->dirichlet_scale * diagu * x[i][j][k].u; in THIFunctionLocal()
824 fn[0]->v = thi->dirichlet_scale * diagv * x[i][j][k].v; in THIFunctionLocal()
827 const PetscReal jw = 0.25 * hx * hy / thi->rhog, *phi = QuadQInterp[q]; in THIFunctionLocal()
835 … THIFriction(thi, PetscRealPart(rbeta2), PetscRealPart(u * u + v * v) / 2, &beta2, &dbeta2); in THIFunctionLocal()
851 PetscCall(PRangeMinMax(&thi->eta, etamin, etamax)); in THIFunctionLocal()
852 PetscCall(PRangeMinMax(&thi->beta2, beta2min, beta2max)); in THIFunctionLocal()
856 static PetscErrorCode THIMatrixStatistics(THI thi, Mat B, PetscViewer viewer) in THIMatrixStatistics() argument
871 … (double)thi->eta.cmin, (double)thi->eta.cmax, (double)thi->beta2.cmin, (double)thi->beta2.cmax)); in THIMatrixStatistics()
904 static PetscErrorCode THISolveStatistics(THI thi, SNES snes, PetscInt coarsened, const char name[]) in THISolveStatistics() argument
911 PetscCall(PetscObjectGetComm((PetscObject)thi, &comm)); in THISolveStatistics()
940 …PetscCallMPI(MPIU_Allreduce(tmin, min, 3, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)thi))); in THISolveStatistics()
941 …PetscCallMPI(MPIU_Allreduce(tmax, max, 3, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)thi))); in THISolveStatistics()
943 nrm2 *= thi->units->year / thi->units->meter; in THISolveStatistics()
945 min[j] *= thi->units->year / thi->units->meter; in THISolveStatistics()
946 max[j] *= thi->units->year / thi->units->meter; in THISolveStatistics()
953 umin *= thi->units->year / thi->units->meter; in THISolveStatistics()
954 umax *= thi->units->year / thi->units->meter; in THISolveStatistics()
955 umean *= thi->units->year / thi->units->meter; in THISolveStatistics()
959 … %g converged range %g to %g\n", (double)thi->eta.min, (double)thi->eta.max, (double)thi->eta.cmin… in THISolveStatistics()
960 … converged range %g to %g\n", (double)thi->beta2.min, (double)thi->beta2.max, (double)thi->beta2.c… in THISolveStatistics()
965 static PetscErrorCode THIJacobianLocal_2D(DMDALocalInfo *info, Node **x, Mat J, Mat B, THI thi) in THIJacobianLocal_2D() argument
976 hx = thi->Lx / info->my; in THIJacobianLocal_2D()
977 hy = thi->Ly / info->mx; in THIJacobianLocal_2D()
1000 jw = 0.25 * hx * hy / thi->rhog; /* rhog is only scaling */ in THIJacobianLocal_2D()
1001 PointwiseNonlinearity2D(thi, n, phi, dphi, &u, &v, du, dv, &eta, &deta); in THIJacobianLocal_2D()
1002 THIFriction(thi, PetscRealPart(rbeta2), PetscRealPart(u * u + v * v) / 2, &beta2, &dbeta2); in THIJacobianLocal_2D()
1011 … * 4. * dpl[0] + dp[1] * jw * eta * dpl[1] + pp * jw * (beta2 / h) * ppl * thi->ssa_friction_scale; in THIJacobianLocal_2D()
1014 … * 4. * dpl[1] + dp[0] * jw * eta * dpl[0] + pp * jw * (beta2 / h) * ppl * thi->ssa_friction_scale; in THIJacobianLocal_2D()
1016 …w * deta * dgdu * (du[1] + dv[0]) + pp * jw * (dbeta2 / h) * u * u * ppl * thi->ssa_friction_scale; in THIJacobianLocal_2D()
1017 …w * deta * dgdv * (du[1] + dv[0]) + pp * jw * (dbeta2 / h) * u * v * ppl * thi->ssa_friction_scale; in THIJacobianLocal_2D()
1018 …w * deta * dgdu * (du[1] + dv[0]) + pp * jw * (dbeta2 / h) * v * u * ppl * thi->ssa_friction_scale; in THIJacobianLocal_2D()
1019 …w * deta * dgdv * (du[1] + dv[0]) + pp * jw * (dbeta2 / h) * v * v * ppl * thi->ssa_friction_scale; in THIJacobianLocal_2D()
1039 if (thi->verbose) PetscCall(THIMatrixStatistics(thi, B, PETSC_VIEWER_STDOUT_WORLD)); in THIJacobianLocal_2D()
1043 static PetscErrorCode THIJacobianLocal_3D(DMDALocalInfo *info, Node ***x, Mat B, THI thi, THIAssemb… argument
1055 hx = thi->Lx / info->mz;
1056 hy = thi->Ly / info->my;
1075 if (thi->no_slip && k == 0) {
1084 PointwiseNonlinearity(thi, n, phi, dphi, &u, &v, du, dv, &eta, &deta);
1085 jw /= thi->rhog; /* residuals are scaled by this factor */
1142 if (thi->no_slip) {
1144 …etscScalar diagu = 2 * etabase / thi->rhog * (hx * hy / hz + hx * hz / hy + 4 * hy * hz / hx), dia…
1145 Ke[0][0] = thi->dirichlet_scale * diagu;
1146 Ke[1][1] = thi->dirichlet_scale * diagv;
1149 const PetscReal jw = 0.25 * hx * hy / thi->rhog, *phi = QuadQInterp[q];
1157 … THIFriction(thi, PetscRealPart(rbeta2), PetscRealPart(u * u + v * v) / 2, &beta2, &dbeta2);
1229 if (thi->verbose) PetscCall(THIMatrixStatistics(thi, B, PETSC_VIEWER_STDOUT_WORLD));
1233 …atic PetscErrorCode THIJacobianLocal_3D_Full(DMDALocalInfo *info, Node ***x, Mat A, Mat B, THI thi) argument
1236 PetscCall(THIJacobianLocal_3D(info, x, B, thi, THIASSEMBLY_FULL));
1240 …tscErrorCode THIJacobianLocal_3D_Tridiagonal(DMDALocalInfo *info, Node ***x, Mat A, Mat B, THI thi) argument
1243 PetscCall(THIJacobianLocal_3D(info, x, B, thi, THIASSEMBLY_TRIDIAGONAL));
1249 THI thi; local
1256 PetscCall(PetscObjectQuery((PetscObject)dac0, "THI", (PetscObject *)&thi));
1257 …PetscCheck(thi, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot refine this DMDA, missing composed T…
1268 …ct)dac), DM_BOUNDARY_NONE, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, st, thi->zlevels, N, M, 1, …
1358 static PetscErrorCode THIDAVecView_VTK_XML(THI thi, DM da, Vec X, const char filename[]) argument
1361 Units units = thi->units;
1369 PetscCall(PetscObjectGetComm((PetscObject)thi, &comm));
1411 PetscReal xx = thi->Lx * i / mx, yy = thi->Ly * j / my, zz;
1412 thi->initialize(thi, xx, yy, &p);
1448 THI thi; local
1456 PetscCall(THICreate(comm, &thi));
1464 if (thi->coarse2d) {
1465 …"-zlevels", "Number of elements in z-direction on fine level", "", thi->zlevels, &thi->zlevels, NU…
1471 if (thi->coarse2d) {
1478 PetscCall(PetscObjectCompose((PetscObject)da, "THI", (PetscObject)thi));
1487 PetscCall(THISetUpDM(thi, da));
1488 if (thi->tridiagonal) da->ops->creatematrix = DMCreateMatrix_THI_Tridiagonal;
1494 if (rlevel - clevel > 0) PetscCall(DMSetMatType(da, thi->mattype));
1497 PetscCall(DMDASNESSetFunctionLocal(da, ADD_VALUES, (DMDASNESFunctionFn *)THIFunctionLocal, thi));
1498 if (thi->tridiagonal) {
1499 …etscCall(DMDASNESSetJacobianLocal(da, (DMDASNESJacobianFn *)THIJacobianLocal_3D_Tridiagonal, thi));
1501 PetscCall(DMDASNESSetJacobianLocal(da, (DMDASNESJacobianFn *)THIJacobianLocal_3D_Full, thi));
1503 PetscCall(DMCoarsenHookAdd(da, DMCoarsenHook_THI, NULL, thi));
1504 PetscCall(DMRefineHookAdd(da, DMRefineHook_THI, NULL, thi));
1506 PetscCall(DMSetApplicationContext(da, thi));
1516 PetscCall(THISolveStatistics(thi, snes, 0, "Full"));
1527 PetscCall(THIDAVecView_VTK_XML(thi, dm, X, filename));
1533 PetscCall(THIDestroy(&thi));