Lines Matching refs:j
285 static inline PetscScalar UInterp(Field **x, PetscInt i, PetscInt j) in UInterp() argument
287 return 0.25 * (x[j][i].u + x[j + 1][i].u + x[j][i + 1].u + x[j + 1][i + 1].u); in UInterp()
290 static inline PetscScalar WInterp(Field **x, PetscInt i, PetscInt j) in WInterp() argument
292 return 0.25 * (x[j][i].w + x[j + 1][i].w + x[j][i + 1].w + x[j + 1][i + 1].w); in WInterp()
295 static inline PetscScalar PInterp(Field **x, PetscInt i, PetscInt j) in PInterp() argument
297 return 0.25 * (x[j][i].p + x[j + 1][i].p + x[j][i + 1].p + x[j + 1][i + 1].p); in PInterp()
300 static inline PetscScalar TInterp(Field **x, PetscInt i, PetscInt j) in TInterp() argument
302 return 0.25 * (x[j][i].T + x[j + 1][i].T + x[j][i + 1].T + x[j + 1][i + 1].T); in TInterp()
306 static inline PetscScalar HorizVelocity(PetscInt i, PetscInt j, AppCtx *user) in HorizVelocity() argument
314 z = (j - grid->jlid - 0.5) * grid->dz; in HorizVelocity()
323 static inline PetscScalar VertVelocity(PetscInt i, PetscInt j, AppCtx *user) in VertVelocity() argument
331 z = (j - grid->jlid) * grid->dz; in VertVelocity()
340 static inline PetscScalar Pressure(PetscInt i, PetscInt j, AppCtx *user) in Pressure() argument
347 z = (j - grid->jlid - 0.5) * grid->dz; in Pressure()
355 static inline PetscScalar CalcSecInv(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user) in CalcSecInv() argument
363 if (i < j) return EPS_ZERO; in CalcSecInv()
365 if (j == jlim) j--; in CalcSecInv()
368 if (j <= grid->jlid) return EPS_ZERO; in CalcSecInv()
370 uE = x[j][i].u; in CalcSecInv()
371 uW = x[j][i - 1].u; in CalcSecInv()
372 wN = x[j][i].w; in CalcSecInv()
373 wS = x[j - 1][i].w; in CalcSecInv()
374 wE = WInterp(x, i, j - 1); in CalcSecInv()
375 if (i == j) { 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()
389 uN = x[j + 1][i].u; in CalcSecInv()
390 uS = x[j][i].u; in CalcSecInv()
391 wE = x[j][i + 1].w; in CalcSecInv()
392 wW = x[j][i].w; in CalcSecInv()
393 if (i == j) { in CalcSecInv()
397 wN = WInterp(x, i, j); in CalcSecInv()
398 uW = UInterp(x, i - 1, j); in CalcSecInv()
401 if (j == grid->jlid) { in CalcSecInv()
407 uE = UInterp(x, i, j); in CalcSecInv()
408 wS = WInterp(x, i, j - 1); in CalcSecInv()
462 static inline PetscScalar XMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user) in XMomentumResidual() argument
475 TS = param->potentialT * TInterp(x, i, j - 1) * PetscExpScalar((j - 1.0) * dz * z_scale); in XMomentumResidual()
476 if (j == jlim) TN = TS; 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()
481 epsN = CalcSecInv(x, i, j, CELL_CORNER, user); in XMomentumResidual()
482 epsS = CalcSecInv(x, i, j - 1, CELL_CORNER, user); in XMomentumResidual()
483 epsE = CalcSecInv(x, i + 1, j, CELL_CENTER, user); in XMomentumResidual()
484 epsW = CalcSecInv(x, i, j, CELL_CENTER, user); in XMomentumResidual()
487 etaN = Viscosity(TN, epsN, dz * (j + 0.5), param); in XMomentumResidual()
488 etaS = Viscosity(TS, epsS, dz * (j - 0.5), param); in XMomentumResidual()
489 etaW = Viscosity(TW, epsW, dz * j, param); in XMomentumResidual()
490 etaE = Viscosity(TE, epsE, dz * j, 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()
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()
513 static inline PetscScalar ZMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user) in ZMomentumResidual() argument
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()
534 epsN = CalcSecInv(x, i, j + 1, CELL_CENTER, user); in ZMomentumResidual()
535 epsS = CalcSecInv(x, i, j, CELL_CENTER, user); in ZMomentumResidual()
536 epsE = CalcSecInv(x, i, j, CELL_CORNER, user); in ZMomentumResidual()
537 epsW = CalcSecInv(x, i - 1, j, CELL_CORNER, user); in ZMomentumResidual()
540 etaN = Viscosity(TN, epsN, dz * (j + 1.0), param); in ZMomentumResidual()
541 etaS = Viscosity(TS, epsS, dz * (j + 0.0), param); in ZMomentumResidual()
542 etaW = Viscosity(TW, epsW, dz * (j + 0.5), param); in ZMomentumResidual()
543 etaE = Viscosity(TE, epsE, dz * (j + 0.5), param); 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()
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()
567 static inline PetscScalar ContinuityResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user) in ContinuityResidual() argument
572 uW = x[j][i - 1].u; in ContinuityResidual()
573 uE = x[j][i].u; in ContinuityResidual()
575 wS = x[j - 1][i].w; in ContinuityResidual()
576 wN = x[j][i].w; in ContinuityResidual()
583 static inline PetscScalar EnergyResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user) in EnergyResidual() argument
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()
602 if (j <= jlid && i >= j) { in EnergyResidual()
605 } else if (i < j) { in EnergyResidual()
611 uW = x[j][i - 1].u; in EnergyResidual()
612 uE = x[j][i].u; in EnergyResidual()
613 wS = x[j - 1][i].w; in EnergyResidual()
614 wN = x[j][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()
620 TN = (x[j][i].T + x[j + 1][i].T) / 2.0; in EnergyResidual()
621 TE = (x[j][i].T + x[j][i + 1].T) / 2.0; in EnergyResidual()
622 TW = (x[j][i].T + x[j][i - 1].T) / 2.0; in EnergyResidual()
630 …[j][i + 2].T + 5.0 * (x[j][i + 1].T + x[j][i].T) - x[j][i - 1].T) / 8.0 - PetscAbsScalar(uE) * (-x… in EnergyResidual()
631 …[j][i + 1].T + 5.0 * (x[j][i].T + x[j][i - 1].T) - x[j][i - 2].T) / 8.0 - PetscAbsScalar(uW) * (-x… in EnergyResidual()
632 …[j + 2][i].T + 5.0 * (x[j + 1][i].T + x[j][i].T) - x[j - 1][i].T) / 8.0 - PetscAbsScalar(wN) * (-x… in EnergyResidual()
633 …[j + 1][i].T + 5.0 * (x[j][i].T + x[j - 1][i].T) - x[j - 2][i].T) / 8.0 - PetscAbsScalar(wS) * (-x… in EnergyResidual()
642 static inline PetscScalar ShearStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *use… in ShearStress() argument
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()
654 if (i == j) { 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()
666 uN = x[j + 1][i].u; in ShearStress()
667 uS = x[j][i].u; in ShearStress()
668 wW = x[j][i].w; in ShearStress()
669 wE = x[j][i + 1].w; in ShearStress()
676 static inline PetscScalar XNormalStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *u… in XNormalStress() argument
683 if (i < j || j <= grid->jlid) return EPS_ZERO; in XNormalStress()
690 TC = param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * z_scale); in XNormalStress()
691 if (ivisc >= VISC_DISL) epsC = CalcSecInv(x, i, j, CELL_CENTER, user); in XNormalStress()
692 etaC = Viscosity(TC, epsC, dz * j, param); in XNormalStress()
694 uW = x[j][i - 1].u; in XNormalStress()
695 uE = x[j][i].u; in XNormalStress()
696 pC = x[j][i].p; in XNormalStress()
699 if (i == ilim || j == jlim) return EPS_ZERO; in XNormalStress()
701 TC = param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * z_scale); in XNormalStress()
702 if (ivisc >= VISC_DISL) epsC = CalcSecInv(x, i, j, CELL_CORNER, user); in XNormalStress()
703 etaC = Viscosity(TC, epsC, dz * (j + 0.5), param); in XNormalStress()
705 if (i == j) uW = param->sb; in XNormalStress()
706 else uW = UInterp(x, i - 1, j); in XNormalStress()
707 uE = UInterp(x, i, j); in XNormalStress()
708 pC = PInterp(x, i, j); in XNormalStress()
715 static inline PetscScalar ZNormalStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *u… in ZNormalStress() argument
723 if (i < j || j <= grid->jlid) return EPS_ZERO; in ZNormalStress()
730 TC = param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * z_scale); in ZNormalStress()
731 if (ivisc >= VISC_DISL) epsC = CalcSecInv(x, i, j, CELL_CENTER, user); in ZNormalStress()
732 etaC = Viscosity(TC, epsC, dz * j, param); in ZNormalStress()
733 wN = x[j][i].w; in ZNormalStress()
734 wS = x[j - 1][i].w; in ZNormalStress()
735 pC = x[j][i].p; in ZNormalStress()
738 if ((i == ilim) || (j == jlim)) return EPS_ZERO; in ZNormalStress()
740 TC = param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * z_scale); in ZNormalStress()
741 if (ivisc >= VISC_DISL) epsC = CalcSecInv(x, i, j, CELL_CORNER, user); in ZNormalStress()
742 etaC = Viscosity(TC, epsC, dz * (j + 0.5), param); in ZNormalStress()
743 if (i == j) wN = param->sb; in ZNormalStress()
744 else wN = WInterp(x, i, j); in ZNormalStress()
745 wS = WInterp(x, i, j - 1); in ZNormalStress()
746 pC = PInterp(x, i, j); in ZNormalStress()
965 PetscInt i, j, is, js, im, jm; in Initialize() local
979 for (j = js; j < js + jm; j++) { 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()
983 else x[j][i].u = HorizVelocity(i, j, user); 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()
987 else x[j][i].w = VertVelocity(i, j, user); in Initialize()
989 if (i < j || j <= grid->jlid) x[j][i].p = 0.0; in Initialize()
990 else x[j][i].p = Pressure(i, j, user); in Initialize()
992 x[j][i].T = PetscMin(grid->dz * (j - 0.5), 1.0); in Initialize()
1100 PetscInt i, j, is, js, im, jm, ilim, jlim, ivt; in ViscosityField() local
1123 for (j = js; j < js + jm; j++) { in ViscosityField()
1125 …T = PetscRealPart(param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * param->z_scale)); in ViscosityField()
1126 if (i < ilim && j < jlim) { in ViscosityField()
1127 …TC = PetscRealPart(param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * param->z_scale)); in ViscosityField()
1131 eps = PetscRealPart(CalcSecInv(x, i, j, CELL_CENTER, user)); in ViscosityField()
1132 epsC = PetscRealPart(CalcSecInv(x, i, j, CELL_CORNER, user)); in ViscosityField()
1134 v[j][i].u = eps; in ViscosityField()
1135 v[j][i].w = epsC; in ViscosityField()
1136 v[j][i].p = Viscosity(T, eps, dz * (j - 0.5), param); in ViscosityField()
1137 v[j][i].T = Viscosity(TC, epsC, dz * j, param); in ViscosityField()
1154 PetscInt i, j, is, js, im, jm; in StressField() local
1171 for (j = js; j < js + jm; j++) { in StressField()
1173 x[j][i].u = ShearStress(y, i, j, CELL_CENTER, user); in StressField()
1174 x[j][i].w = ShearStress(y, i, j, CELL_CORNER, user); in StressField()
1175 x[j][i].p = XNormalStress(y, i, j, CELL_CENTER, user); in StressField()
1176 x[j][i].T = ZNormalStress(y, i, j, CELL_CENTER, user); in StressField()
1192 static inline PetscScalar SlabVel(char c, PetscInt i, PetscInt j, AppCtx *user) in SlabVel() argument
1198 if (i < j - 1) return param->cb; in SlabVel()
1199 else if (j <= grid->jfault) return 0.0; in SlabVel()
1203 if (i < j) return param->sb; in SlabVel()
1204 else if (j <= grid->jfault) return 0.0; in SlabVel()
1210 static inline PetscScalar PlateModel(PetscInt j, PetscInt plate, AppCtx *user) in PlateModel() argument
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()
1294 PetscInt i, j, mx, mz, ilim, jlim; in FormFunctionLocal() local
1311 for (j = js; j < je; j++) { 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()
1326 } else if (j == jlim) { in FormFunctionLocal()
1329 f[j][i].u = x[j][i].u - HorizVelocity(i, j, user); in FormFunctionLocal()
1331 f[j][i].u = XMomentumResidual(x, i, j, user); in FormFunctionLocal()
1338 f[j][i].u = XMomentumResidual(x, i, j, user); in FormFunctionLocal()
1342 if (i <= j) { 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()
1349 } else if (j == jlim) { 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()
1362 f[j][i].w = ZMomentumResidual(x, i, j, user); in FormFunctionLocal()
1369 f[j][i].w = ZMomentumResidual(x, i, j, user); in FormFunctionLocal()
1373 …if (i < j || j <= grid->jlid || (j < grid->corner + grid->inose && i < grid->corner + grid->inose)… in FormFunctionLocal()
1375 f[j][i].p = x[j][i].p; in FormFunctionLocal()
1377 } else if ((i == ilim || j == jlim) && ibound == BC_ANALYTIC) { in FormFunctionLocal()
1379 f[j][i].p = x[j][i].p - Pressure(i, j, user); in FormFunctionLocal()
1383 f[j][i].p = ContinuityResidual(x, i, j, user); in FormFunctionLocal()
1387 if (j == 0) { in FormFunctionLocal()
1389 f[j][i].T = x[j][i].T + x[j + 1][i].T + PetscMax(PetscRealPart(x[j][i].T), 0.0); 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()
1400 } else if (j == jlim) { 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()
1407 f[j][i].T = EnergyResidual(x, i, j, user); in FormFunctionLocal()