435 xt::pyarray<double> &mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
436 xt::pyarray<double> &mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
437 xt::pyarray<double> &mesh_dof = args.
array<
double>(
"mesh_dof");
438 xt::pyarray<double> &mesh_velocity_dof = args.
array<
double>(
"mesh_velocity_dof");
439 double MOVING_DOMAIN = args.
scalar<
double>(
"MOVING_DOMAIN");
440 xt::pyarray<int> &mesh_l2g = args.
array<
int>(
"mesh_l2g");
441 xt::pyarray<double> &dV_ref = args.
array<
double>(
"dV_ref");
442 xt::pyarray<double> &u_trial_ref = args.
array<
double>(
"u_trial_ref");
443 xt::pyarray<double> &u_grad_trial_ref = args.
array<
double>(
"u_grad_trial_ref");
444 xt::pyarray<double> &u_test_ref = args.
array<
double>(
"u_test_ref");
445 xt::pyarray<double> &u_grad_test_ref = args.
array<
double>(
"u_grad_test_ref");
446 xt::pyarray<double> &mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
447 xt::pyarray<double> &mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
448 xt::pyarray<double> &dS_ref = args.
array<
double>(
"dS_ref");
449 xt::pyarray<double> &u_trial_trace_ref = args.
array<
double>(
"u_trial_trace_ref");
450 xt::pyarray<double> &u_grad_trial_trace_ref = args.
array<
double>(
"u_grad_trial_trace_ref");
451 xt::pyarray<double> &u_test_trace_ref = args.
array<
double>(
"u_test_trace_ref");
452 xt::pyarray<double> &u_grad_test_trace_ref = args.
array<
double>(
"u_grad_test_trace_ref");
453 xt::pyarray<double> &normal_ref = args.
array<
double>(
"normal_ref");
454 xt::pyarray<double> &boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
455 int nElements_global = args.
scalar<
int>(
"nElements_global");
456 xt::pyarray<double> &ebqe_penalty_ext = args.
array<
double>(
"ebqe_penalty_ext");
457 xt::pyarray<int> &elementMaterialTypes = args.
array<
int>(
"elementMaterialTypes");
458 xt::pyarray<int> &isSeepageFace = args.
array<
int>(
"isSeepageFace");
459 xt::pyarray<int> &a_rowptr = args.
array<
int>(
"a_rowptr");
460 xt::pyarray<int> &a_colind = args.
array<
int>(
"a_colind");
461 double rho = args.
scalar<
double>(
"rho");
462 double beta = args.
scalar<
double>(
"beta");
465 xt::pyarray<double> &q_rho = args.
array<
double>(
"q_rho");
466 xt::pyarray<double> &ebqe_rho = args.
array<
double>(
"ebqe_rho");
468 xt::pyarray<double> &gravity = args.
array<
double>(
"gravity");
469 xt::pyarray<double> &alpha = args.
array<
double>(
"alpha");
470 xt::pyarray<double> &
n = args.
array<
double>(
"n");
471 xt::pyarray<double> &thetaR = args.
array<
double>(
"thetaR");
472 xt::pyarray<double> &thetaSR = args.
array<
double>(
"thetaSR");
473 xt::pyarray<double> &KWs = args.
array<
double>(
"KWs");
474 xt::pyarray<double> &krn_end = args.
array<
double>(
"krn_end");
475 xt::pyarray<double> &S_gr = args.
array<
double>(
"S_gr");
476 double mu_n = args.
scalar<
double>(
"mu_n");
477 double useMetrics = args.
scalar<
double>(
"useMetrics");
478 double alphaBDF = args.
scalar<
double>(
"alphaBDF");
479 int lag_shockCapturing = args.
scalar<
int>(
"lag_shockCapturing");
480 double shockCapturingDiffusion = args.
scalar<
double>(
"shockCapturingDiffusion");
481 double sc_uref = args.
scalar<
double>(
"sc_uref");
482 double sc_alpha = args.
scalar<
double>(
"sc_alpha");
483 xt::pyarray<int> &u_l2g = args.
array<
int>(
"u_l2g");
484 xt::pyarray<double> &elementDiameter = args.
array<
double>(
"elementDiameter");
485 xt::pyarray<double> &u_dof = args.
array<
double>(
"u_dof");
486 xt::pyarray<double> &u_dof_old = args.
array<
double>(
"u_dof_old");
487 xt::pyarray<double> &velocity = args.
array<
double>(
"velocity");
488 xt::pyarray<double> &q_m = args.
array<
double>(
"q_m");
489 xt::pyarray<double> &q_theta = args.
array<
double>(
"q_theta");
490 xt::pyarray<double> &q_u = args.
array<
double>(
"q_u");
491 xt::pyarray<double> &q_dV = args.
array<
double>(
"q_dV");
492 xt::pyarray<double> &q_m_betaBDF = args.
array<
double>(
"q_m_betaBDF");
493 xt::pyarray<double> &cfl = args.
array<
double>(
"cfl");
494 xt::pyarray<double> &q_numDiff_u = args.
array<
double>(
"q_numDiff_u");
495 xt::pyarray<double> &q_numDiff_u_last = args.
array<
double>(
"q_numDiff_u_last");
496 int offset_u = args.
scalar<
int>(
"offset_u");
497 int stride_u = args.
scalar<
int>(
"stride_u");
501 const double dt = args.
scalar<
double>(
"dt");
502 xt::pyarray<double> &u_dof_n = args.
array<
double>(
"u_dof_n");
503 xt::pyarray<double> &u_dof_n_old = args.
array<
double>(
"u_dof_n_old");
507 const double rho_n = args.
scalar<
double>(
"rho_n");
508 const double p_ref_n = args.
scalar<
double>(
"p_ref_n");
509 const bool rho_n_compressible = (p_ref_n > 0.0);
510 const double c_n = rho_n_compressible ? (rho_n / p_ref_n) : 0.0;
511 const int offset_n = args.
scalar<
int>(
"offset_n");
512 const int stride_n = args.
scalar<
int>(
"stride_n");
518 xt::pyarray<double> &c_dof = args.
array<
double>(
"c_dof");
519 const double k_d = args.
scalar<
double>(
"k_d");
520 const double c_sat = args.
scalar<
double>(
"c_sat");
524 xt::pyarray<double> &injection_dof = args.
array<
double>(
"injection_dof");
525 xt::pyarray<double> &globalResidual = args.
array<
double>(
"globalResidual");
526 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
527 xt::pyarray<int> &exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
528 xt::pyarray<int> &elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
529 xt::pyarray<int> &elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
530 xt::pyarray<double> &ebqe_velocity_ext = args.
array<
double>(
"ebqe_velocity_ext");
531 xt::pyarray<int> &isDOFBoundary_u = args.
array<
int>(
"isDOFBoundary_u");
532 xt::pyarray<double> &ebqe_bc_u_ext = args.
array<
double>(
"ebqe_bc_u_ext");
534 xt::pyarray<int> &isDOFBoundary_n = args.
array<
int>(
"isDOFBoundary_n");
535 xt::pyarray<double> &ebqe_bc_u_n_ext = args.
array<
double>(
"ebqe_bc_u_n_ext");
536 xt::pyarray<int> &isFluxBoundary_u = args.
array<
int>(
"isFluxBoundary_u");
537 xt::pyarray<double> &ebqe_bc_flux_ext = args.
array<
double>(
"ebqe_bc_flux_ext");
538 xt::pyarray<double> &ebqe_phi = args.
array<
double>(
"ebqe_phi");
539 double epsFact = args.
scalar<
double>(
"epsFact");
540 xt::pyarray<double> &ebqe_u = args.
array<
double>(
"ebqe_u");
541 xt::pyarray<double> &ebqe_theta = args.
array<
double>(
"ebqe_theta");
542 xt::pyarray<double> &ebqe_flux = args.
array<
double>(
"ebqe_flux");
546 double cE = args.
scalar<
double>(
"cE");
547 double cK = args.
scalar<
double>(
"cK");
549 double uL = args.
scalar<
double>(
"uL");
550 double uR = args.
scalar<
double>(
"uR");
552 int numDOFs = args.
scalar<
int>(
"numDOFs");
556 int numDOFs_u = args.
scalar<
int>(
"numDOFs_u");
557 int NNZ = args.
scalar<
int>(
"NNZ");
558 xt::pyarray<int> &csrRowIndeces_DofLoops = args.
array<
int>(
"csrRowIndeces_DofLoops");
559 xt::pyarray<int> &csrColumnOffsets_DofLoops = args.
array<
int>(
"csrColumnOffsets_DofLoops");
560 xt::pyarray<int> &csrRowIndeces_Full = args.
array<
int>(
"csrRowIndeces_Full");
561 xt::pyarray<int> &csrColumnOffsets_Full = args.
array<
int>(
"csrColumnOffsets_Full");
562 xt::pyarray<int> &csrRowIndeces_CellLoops = args.
array<
int>(
"csrRowIndeces_CellLoops");
563 xt::pyarray<int> &csrColumnOffsets_CellLoops = args.
array<
int>(
"csrColumnOffsets_CellLoops");
564 xt::pyarray<int> &csrColumnOffsets_eb_CellLoops = args.
array<
int>(
"csrColumnOffsets_eb_CellLoops");
566 xt::pyarray<double> &Cx = args.
array<
double>(
"Cx");
567 xt::pyarray<double> &Cy = args.
array<
double>(
"Cy");
568 xt::pyarray<double> &Cz = args.
array<
double>(
"Cz");
569 xt::pyarray<double> &CTx = args.
array<
double>(
"CTx");
570 xt::pyarray<double> &CTy = args.
array<
double>(
"CTy");
571 xt::pyarray<double> &CTz = args.
array<
double>(
"CTz");
572 xt::pyarray<double> &ML = args.
array<
double>(
"ML");
573 xt::pyarray<double> &delta_x_ij = args.
array<
double>(
"delta_x_ij");
575 int LUMPED_MASS_MATRIX = args.
scalar<
int>(
"LUMPED_MASS_MATRIX");
577 int ENTROPY_TYPE = args.
scalar<
int>(
"ENTROPY_TYPE");
583 xt::pyarray<double> &dLow = args.
array<
double>(
"dLow");
584 xt::pyarray<double> &fluxMatrix = args.
array<
double>(
"fluxMatrix");
586 xt::pyarray<double> &quantDOFs = args.
array<
double>(
"quantDOFs");
588 assert(a_rowptr.data()[nSpace] ==
nnz);
589 assert(a_rowptr.data()[nSpace] == nSpace);
593 xt::pyarray<double> &anb_seepage_flux_n = args.
array<
double>(
"anb_seepage_flux_n");
595 xt::pyarray<double> &velocity_couple = args.
array<
double>(
"velocity_couple");
596 xt::pyarray<double> &ebqe_velocity_ext_couple = args.
array<
double>(
"ebqe_velocity_ext_couple");
602 double &anb_seepage_flux(args.
scalar<
double>(
"anb_seepage_flux"));
603 xt::pyarray<double> &q_velocity = args.
array<
double>(
"q_velocity");
604 anb_seepage_flux = 0.0;
615 for (
int eN = 0; eN < nElements_global; eN++) {
617 double elementResidual_u[nDOF_test_element];
618 for (
int i = 0; i < nDOF_test_element; i++) { elementResidual_u[i] = 0.0; }
620 for (
int k = 0; k < nQuadraturePoints_element; k++) {
622 int eN_k = eN * nQuadraturePoints_element + k, eN_k_nSpace = eN_k * nSpace, eN_nDOF_trial_element = eN * nDOF_trial_element;
623 double u = 0.0, grad_u[nSpace], grad_u_old[nSpace], m = 0.0, dm = 0.0,
f[nSpace],
df[nSpace], a[
nnz], da[
nnz], as[
nnz], m_t = 0.0, dm_t = 0.0, pdeResidual_u = 0.0, Lstar_u[nDOF_test_element], subgridError_u = 0.0, tau = 0.0, tau0 = 0.0, tau1 = 0.0, numDiff0 = 0.0, numDiff1 = 0.0, jac[nSpace * nSpace], jacDet, jacInv[nSpace * nSpace], u_grad_trial[nDOF_trial_element * nSpace], u_test_dV[nDOF_trial_element], u_grad_test_dV[nDOF_test_element * nSpace], dV, x, y,
z,
xt, yt, zt, G[nSpace * nSpace], G_dd_G, tr_G, norm_Rv;
627 ck.calculateMapping_element(eN, k, mesh_dof.data(), mesh_l2g.data(), mesh_trial_ref.data(), mesh_grad_trial_ref.data(), jac, jacDet, jacInv, x, y,
z);
628 ck.calculateMappingVelocity_element(eN, k, mesh_velocity_dof.data(), mesh_l2g.data(), mesh_trial_ref.data(),
xt, yt, zt);
630 dV = fabs(jacDet) * dV_ref.data()[k];
631 q_dV.data()[eN_k] = dV;
632 ck.calculateG(jacInv, G, G_dd_G, tr_G);
634 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k * nDOF_trial_element * nSpace], jacInv, u_grad_trial);
636 ck.valFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], &u_trial_ref.data()[k * nDOF_trial_element],
u);
638 ck.gradFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], u_grad_trial, grad_u);
641 double grad_u_n[nSpace];
642 ck.gradFromDOF(u_dof_n.data(), &u_l2g.data()[eN_nDOF_trial_element], u_grad_trial, grad_u_n);
651 for (
int j = 0; j < nDOF_trial_element; j++) {
652 u_test_dV[j] = u_test_ref.data()[k * nDOF_trial_element + j] * dV;
653 for (
int I = 0; I < nSpace; I++) {
654 u_grad_test_dV[j * nSpace + I] = u_grad_trial[j * nSpace + I] * dV;
660 double Kr, dKr, thetaW;
661 const double rho_local = q_rho.data()[eN_k];
662 const double rho_velocity = std::fabs(rho_local) > 1.0e-12 ? rho_local : rho;
665 double dm_du_n_qp = 0.0, dkr_du_n_qp = 0.0;
666 double df_du_n_qp[nSpace];
667 double da_du_n_qp[
nnz];
668 for (
int I = 0; I < nSpace; I++) df_du_n_qp[I] = 0.0;
669 for (
int ii = 0; ii <
nnz; ii++) da_du_n_qp[ii] = 0.0;
671 ck.valFromDOF(u_dof_n.data(),
672 &u_l2g.data()[eN_nDOF_trial_element],
673 &u_trial_ref.data()[k * nDOF_trial_element], u_n_qp);
675 alpha.data()[elementMaterialTypes.data()[eN]],
n.data()[elementMaterialTypes.data()[eN]],
676 thetaR.data()[elementMaterialTypes.data()[eN]], thetaSR.data()[elementMaterialTypes.data()[eN]],
677 &KWs.data()[elementMaterialTypes.data()[eN] *
nnz],
u, u_n_qp,
678 m, dm, dm_du_n_qp,
f,
df, df_du_n_qp, a, da, da_du_n_qp,
679 as, Kr, dKr, dkr_du_n_qp, thetaW);
680 q_theta.data()[eN_k] = thetaW;
683 for (
int I = 0; I < nSpace; ++I) {
684 q_velocity.data()[eN_k_nSpace + I] = grad_u[I];
687 double pressure_gradient[nSpace];
688 for (
int J=0; J<nSpace; ++J)
689 pressure_gradient[J] = grad_u[J] - rho_velocity * gravity.data()[J];
691 for (
int I=0; I<nSpace; ++I) {
693 for (
int ii = a_rowptr.data()[I]; ii < a_rowptr.data()[I+1]; ++ii) {
694 const int J = a_colind.data()[ii];
695 acc += (a[ii] / rho_velocity) * pressure_gradient[J];
697 velocity.data()[eN_k_nSpace + I] = -acc;
698 velocity_couple.data()[eN_k_nSpace + I] = -acc ;
705 const int mat0 = elementMaterialTypes.data()[eN];
706 const double phi0_m = thetaR.data()[mat0] + thetaSR.data()[mat0];
707 const double z0_m = fmin(fmax(u_n_qp, 1.0e-8), 1.0 - 1.0e-8);
708 const double p0_m = fmax(
u, 1.0e2);
714 m = phi0_m * Nm * (1.0 - z0_m);
715 dm = phi0_m * dNm_dp * (1.0 - z0_m);
720 ck.bdf(alphaBDF, q_m_betaBDF.data()[eN_k], m, dm, m_t, dm_t);
748 const int mat_eN0 = elementMaterialTypes.data()[eN];
749 const double alpha_eN0 = alpha.data()[mat_eN0];
750 const double n_vg_eN0 =
n.data()[mat_eN0];
751 const double krn_end0 = krn_end.data()[mat_eN0];
752 const double *KWs_eN0 = &KWs.data()[mat_eN0 *
nnz];
753 const double phi0 = thetaR.data()[mat_eN0] + thetaSR.data()[mat_eN0];
754 const double S_wr0 = thetaR.data()[mat_eN0] / phi0;
755 const double one_m_Sr0 = 1.0 - S_wr0;
756 const double Se_trap_L771 = 1.0 - S_gr.data()[mat_eN0] / one_m_Sr0;
757 const double z_cl0 = fmin(fmax(u_n_qp, 1.0e-8), 1.0 - 1.0e-8);
758 const double p_cl0 = fmax(
u, 1.0e2);
761 const double Se_a0 = fmin(fmax((1.0 - fs0.
S_g - S_wr0) / one_m_Sr0, 0.0), 1.0);
762 double KWr0 = 0.0, DKWr0 = 0.0, thW0 = 0.0, DthW0 = 0.0;
763 double KNr0 = 0.0, DKNr0 = 0.0, pc0 = 0.0, dpc_dSe0 = 0.0, d2pc0 = 0.0;
766 thetaR.data()[mat_eN0], thetaSR.data()[mat_eN0], thW0, DthW0, KWr0, DKWr0);
771 thetaR.data()[mat_eN0], thetaSR.data()[mat_eN0], thW0, DthW0, KWr0, DKWr0);
776 const double pcp0 = dpc_dSe0 / one_m_Sr0;
779 const double rho_g_mass0 = fs0.
rho_g*Mbar_g0;
780 const double rho_a_mass0 = fs0.
rho_a*Mbar_a0;
782 for (
int I = 0; I < nSpace; I++) {
783 double ua = 0.0, ug = 0.0;
784 for (
int ii = a_rowptr.data()[I]; ii < a_rowptr.data()[I + 1]; ii++) {
785 const int J = a_colind.data()[ii];
786 const double gradSa_J = -(fs0.
dS_g_dp*grad_u[J] + fs0.
dS_g_dz*grad_u_n[J]);
787 const double gp_a = grad_u[J] - rho_a_mass0*gravity.data()[J];
788 const double gp_g = grad_u[J] + pcp0*gradSa_J - rho_g_mass0*gravity.data()[J];
789 ua -= (KWr0*KWs_eN0[ii]) * gp_a;
790 ug -= (KNr0*KWs_eN0[ii]/mu_n) * gp_g;
792 F0[I] = fs0.
rho_g*(1.0 - fs0.
Y)*ug + fs0.
rho_a*(1.0 - fs0.
X)*ua;
797 for (
int i = 0; i < nDOF_test_element; i++) {
798 int eN_k_i = eN_k * nDOF_test_element + i, eN_k_i_nSpace = eN_k_i * nSpace, i_nSpace = i * nSpace;
800 elementResidual_u[i] +=
ck.Mass_weak(m_t, u_test_dV[i]) +
VMS *
ck.SubgridError(subgridError_u, Lstar_u[i]) +
VMS *
ck.NumericalDiffusion(q_numDiff_u_last[eN_k], grad_u, &u_grad_test_dV[i_nSpace]);
801 for (
int I = 0; I < nSpace; I++)
802 elementResidual_u[i] -= F0[I] * u_grad_test_dV[i_nSpace + I];
805 q_m.data()[eN_k] = m;
806 q_u.data()[eN_k] =
u;
811 for (
int i = 0; i < nDOF_test_element; i++) {
812 int eN_i = eN * nDOF_test_element + i;
814 globalResidual.data()[offset_u + stride_u * u_l2g.data()[eN_i]] += elementResidual_u[i];
823 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++) {
824 int ebN = exteriorElementBoundariesArray.data()[ebNE], eN = elementBoundaryElementsArray.data()[ebN * 2 + 0], ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN * 2 + 0], eN_nDOF_trial_element = eN * nDOF_trial_element;
825 double elementResidual_u[nDOF_test_element];
826 for (
int i = 0; i < nDOF_test_element; i++) { elementResidual_u[i] = 0.0; }
827 for (
int kb = 0; kb < nQuadraturePoints_elementBoundary; kb++) {
828 int ebNE_kb = ebNE * nQuadraturePoints_elementBoundary + kb, ebNE_kb_nSpace = ebNE_kb * nSpace, ebN_local_kb = ebN_local * nQuadraturePoints_elementBoundary + kb, ebN_local_kb_nSpace = ebN_local_kb * nSpace;
829 double u_ext = 0.0, grad_u_ext[nSpace], m_ext = 0.0, dm_ext = 0.0, f_ext[nSpace], df_ext[nSpace], a_ext[
nnz], da_ext[
nnz], as_ext[
nnz], flux_ext = 0.0,
831 bc_u_ext = 0.0, bc_grad_u_ext[nSpace], bc_m_ext = 0.0, bc_dm_ext = 0.0, bc_f_ext[nSpace], bc_df_ext[nSpace], bc_a_ext[
nnz], bc_da_ext[
nnz], bc_as_ext[
nnz], jac_ext[nSpace * nSpace], jacDet_ext, jacInv_ext[nSpace * nSpace], boundaryJac[nSpace * (nSpace - 1)], metricTensor[(nSpace - 1) * (nSpace - 1)], metricTensorDetSqrt, dS, u_test_dS[nDOF_test_element], u_grad_trial_trace[nDOF_trial_element * nSpace], normal[3], x_ext, y_ext, z_ext, xt_ext, yt_ext, zt_ext, integralScaling, G[nSpace * nSpace], G_dd_G, tr_G;
836 ck.calculateMapping_elementBoundary(eN, ebN_local, kb, ebN_local_kb, mesh_dof.data(), mesh_l2g.data(), mesh_trial_trace_ref.data(), mesh_grad_trial_trace_ref.data(), boundaryJac_ref.data(), jac_ext, jacDet_ext, jacInv_ext, boundaryJac, metricTensor, metricTensorDetSqrt,
837 normal_ref.data(), normal, x_ext, y_ext, z_ext);
838 ck.calculateMappingVelocity_elementBoundary(eN, ebN_local, kb, ebN_local_kb, mesh_velocity_dof.data(), mesh_l2g.data(), mesh_trial_trace_ref.data(), xt_ext, yt_ext, zt_ext, normal, boundaryJac, metricTensor, integralScaling);
839 dS = ((1.0 - MOVING_DOMAIN) * metricTensorDetSqrt + MOVING_DOMAIN * integralScaling) * dS_ref.data()[kb];
842 ck.calculateG(jacInv_ext, G, G_dd_G, tr_G);
845 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace * nDOF_trial_element], jacInv_ext, u_grad_trial_trace);
847 ck.valFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], &u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element], u_ext);
848 ck.gradFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], u_grad_trial_trace, grad_u_ext);
857 for (
int j = 0; j < nDOF_trial_element; j++) { u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb * nDOF_test_element + j] * dS; }
861 bc_u_ext = isDOFBoundary_u.data()[ebNE_kb] * ebqe_bc_u_ext.data()[ebNE_kb] + (1 - isDOFBoundary_u.data()[ebNE_kb]) * u_ext;
865 const double rho_ext = ebqe_rho.data()[ebNE_kb];
866 const double rho_velocity_ext = std::fabs(rho_ext) > 1.0e-12 ? rho_ext : rho;
867 double Kr, dKr, thetaW_ext, thetaW_bc;
870 double dm_du_n_ext = 0.0, dkr_du_n_ext = 0.0;
871 double df_du_n_ext[nSpace];
872 double da_du_n_ext[
nnz];
873 double bc_dm_du_n = 0.0, bc_dkr_du_n = 0.0;
874 double bc_df_du_n[nSpace];
875 double bc_da_du_n[
nnz];
876 for (
int I = 0; I < nSpace; I++) { df_du_n_ext[I] = 0.0; bc_df_du_n[I] = 0.0; }
877 for (
int ii = 0; ii <
nnz; ii++) { da_du_n_ext[ii] = 0.0; bc_da_du_n[ii] = 0.0; }
878 double u_n_ext_qp = 0.0;
879 ck.valFromDOF(u_dof_n.data(), &u_l2g.data()[eN_nDOF_trial_element],
880 &u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element], u_n_ext_qp);
883 const double bc_u_n_ext_qp = isDOFBoundary_n.data()[ebNE_kb] * ebqe_bc_u_n_ext.data()[ebNE_kb]
884 + (1 - isDOFBoundary_n.data()[ebNE_kb]) * u_n_ext_qp;
886 alpha.data()[elementMaterialTypes.data()[eN]],
n.data()[elementMaterialTypes.data()[eN]],
887 thetaR.data()[elementMaterialTypes.data()[eN]], thetaSR.data()[elementMaterialTypes.data()[eN]],
888 &KWs.data()[elementMaterialTypes.data()[eN] *
nnz], u_ext, u_n_ext_qp,
889 m_ext, dm_ext, dm_du_n_ext, f_ext, df_ext, df_du_n_ext, a_ext, da_ext, da_du_n_ext,
890 as_ext, Kr, dKr, dkr_du_n_ext, thetaW_ext);
892 alpha.data()[elementMaterialTypes.data()[eN]],
n.data()[elementMaterialTypes.data()[eN]],
893 thetaR.data()[elementMaterialTypes.data()[eN]], thetaSR.data()[elementMaterialTypes.data()[eN]],
894 &KWs.data()[elementMaterialTypes.data()[eN] *
nnz], bc_u_ext, bc_u_n_ext_qp,
895 bc_m_ext, bc_dm_ext, bc_dm_du_n, bc_f_ext, bc_df_ext, bc_df_du_n, bc_a_ext, bc_da_ext, bc_da_du_n,
896 bc_as_ext, Kr, dKr, bc_dkr_du_n, thetaW_bc);
897 ebqe_theta.data()[ebNE_kb] = thetaW_ext;
902 double ext_pressure_gradient[nSpace];
903 for (
int J=0; J<nSpace; ++J)
904 ext_pressure_gradient[J] = grad_u_ext[J] - rho_velocity_ext * gravity.data()[J];
906 for (
int I=0; I<nSpace; ++I) {
908 for (
int ii = a_rowptr.data()[I]; ii < a_rowptr.data()[I+1]; ++ii) {
909 const int J = a_colind.data()[ii];
910 acc += (a_ext[ii] / rho_velocity_ext) * ext_pressure_gradient[J];
912 ebqe_velocity_ext.data()[ebNE_kb_nSpace + I] = -acc;
913 ebqe_velocity_ext_couple.data()[ebNE_kb_nSpace + I] = -acc ;
923 double grad_u_n_ext[nSpace];
924 ck.gradFromDOF(u_dof_n.data(), &u_l2g.data()[eN_nDOF_trial_element],
925 u_grad_trial_trace, grad_u_n_ext);
926 const int mat_b = elementMaterialTypes.data()[eN];
927 const double alpha_b = alpha.data()[mat_b];
928 const double n_vg_b =
n.data()[mat_b];
929 const double krn_end_b = krn_end.data()[mat_b];
930 const double *KWs_b = &KWs.data()[mat_b *
nnz];
931 const double phi_b = thetaR.data()[mat_b] + thetaSR.data()[mat_b];
932 const double S_wr_b = thetaR.data()[mat_b] / phi_b;
933 const double one_m_Sr_b = 1.0 - S_wr_b;
934 const double Se_trap_L956 = 1.0 - S_gr.data()[mat_b] / one_m_Sr_b;
935 const double z_clb = fmin(fmax(u_n_ext_qp, 1.0e-8), 1.0 - 1.0e-8);
936 const double p_clb = fmax(u_ext, 1.0e2);
939 const double Se_ab = fmin(fmax((1.0 - fsb.
S_g - S_wr_b)/one_m_Sr_b, 0.0), 1.0);
940 double KWrb=0,DKWrb=0,thWb=0,DthWb=0,KNrb=0,DKNrb=0,pcb=0,dpc_dSeb=0,d2pcb=0;
951 const double pcpb = dpc_dSeb / one_m_Sr_b;
954 const double rho_g_mass_b = fsb.
rho_g*Mbar_gb;
955 const double rho_a_mass_b = fsb.
rho_a*Mbar_ab;
957 for (
int I = 0; I < nSpace; I++) {
958 double ua = 0.0, ug = 0.0;
959 for (
int ii = a_rowptr.data()[I]; ii < a_rowptr.data()[I+1]; ii++) {
960 const int J = a_colind.data()[ii];
961 const double gradSa_J = -(fsb.
dS_g_dp*grad_u_ext[J] + fsb.
dS_g_dz*grad_u_n_ext[J]);
962 const double gp_a = grad_u_ext[J] - rho_a_mass_b*gravity.data()[J];
963 const double gp_g = grad_u_ext[J] + pcpb*gradSa_J - rho_g_mass_b*gravity.data()[J];
964 ua -= (KWrb*KWs_b[ii]) * gp_a;
965 ug -= (KNrb*KWs_b[ii]/mu_n) * gp_g;
967 F0n += (fsb.
rho_g*(1.0 - fsb.
Y)*ug + fsb.
rho_a*(1.0 - fsb.
X)*ua) * normal[I];
969 const int isSeep = isSeepageFace.data()[ebNE];
970 if (isSeep || isDOFBoundary_u.data()[ebNE_kb]) {
971 const double bc_u_pen = isSeep ? 0.0 : bc_u_ext;
972 flux_ext = F0n + ebqe_penalty_ext.data()[ebNE_kb] * (u_ext - bc_u_pen);
973 if (isSeep && flux_ext <= 0.0) flux_ext = 0.0;
975 flux_ext = ebqe_bc_flux_ext[ebNE_kb];
978 ebqe_flux.data()[ebNE_kb] = flux_ext;
980 anb_seepage_flux =
seepagefluxcalculator(anb_seepage_flux, isSeepageFace.data()[ebNE], dS, flux_ext);
981 anb_seepage_flux_n.data()[0] = anb_seepage_flux;
982 ebqe_u.data()[ebNE_kb] = u_ext;
986 for (
int i = 0; i < nDOF_test_element; i++) {
987 elementResidual_u[i] +=
ck.ExteriorElementBoundaryFlux(flux_ext, u_test_dS[i]);
994 for (
int i = 0; i < nDOF_test_element; i++) {
995 int eN_i = eN * nDOF_test_element + i;
996 globalResidual.data()[offset_u + stride_u * u_l2g.data()[eN_i]] += elementResidual_u[i];
1016 for (
int eN = 0; eN < nElements_global; eN++) {
1017 const int mat_eN = elementMaterialTypes.data()[eN];
1018 const double phi_eN = thetaR.data()[mat_eN] + thetaSR.data()[mat_eN];
1019 const double alpha_eN = alpha.data()[mat_eN];
1020 const double krn_end_eN = krn_end.data()[mat_eN];
1021 const double n_vg_eN =
n.data()[mat_eN];
1022 const double *KWs_eN = &KWs.data()[mat_eN *
nnz];
1023 double elementResidual_n[nDOF_test_element];
1029 double lumped_w_n[nDOF_test_element];
1030 for (
int i = 0; i < nDOF_test_element; i++) {
1031 elementResidual_n[i] = 0.0;
1032 lumped_w_n[i] = 0.0;
1034 for (
int k = 0; k < nQuadraturePoints_element; k++) {
1035 const int eN_k = eN * nQuadraturePoints_element + k;
1036 const int eN_nDOF_trial_element = eN * nDOF_trial_element;
1037 double jac[nSpace * nSpace], jacDet, jacInv[nSpace * nSpace], x_q, y_q, z_q;
1038 ck.calculateMapping_element(eN, k, mesh_dof.data(), mesh_l2g.data(),
1039 mesh_trial_ref.data(), mesh_grad_trial_ref.data(),
1040 jac, jacDet, jacInv, x_q, y_q, z_q);
1041 const double dV = std::fabs(jacDet) * dV_ref.data()[k];
1043 double u_grad_trial_qp[nDOF_trial_element * nSpace];
1044 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k * nDOF_trial_element * nSpace],
1045 jacInv, u_grad_trial_qp);
1047 double u_n = 0.0, u_n_old = 0.0;
1048 ck.valFromDOF(u_dof_n.data(),
1049 &u_l2g.data()[eN_nDOF_trial_element],
1050 &u_trial_ref.data()[k * nDOF_trial_element], u_n);
1051 ck.valFromDOF(u_dof_n_old.data(),
1052 &u_l2g.data()[eN_nDOF_trial_element],
1053 &u_trial_ref.data()[k * nDOF_trial_element], u_n_old);
1056 double u_w_qp = 0.0, u_w_qp_old = 0.0;
1057 ck.valFromDOF(u_dof.data(),
1058 &u_l2g.data()[eN_nDOF_trial_element],
1059 &u_trial_ref.data()[k * nDOF_trial_element], u_w_qp);
1060 ck.valFromDOF(u_dof_old.data(),
1061 &u_l2g.data()[eN_nDOF_trial_element],
1062 &u_trial_ref.data()[k * nDOF_trial_element], u_w_qp_old);
1063 double grad_u_w[nSpace];
1064 ck.gradFromDOF(u_dof.data(),
1065 &u_l2g.data()[eN_nDOF_trial_element],
1066 u_grad_trial_qp, grad_u_w);
1068 double grad_u_n[nSpace];
1069 ck.gradFromDOF(u_dof_n.data(),
1070 &u_l2g.data()[eN_nDOF_trial_element],
1071 u_grad_trial_qp, grad_u_n);
1079 const double phi_loc = thetaR.data()[mat_eN] + thetaSR.data()[mat_eN];
1080 const double S_wr_loc = thetaR.data()[mat_eN] / phi_loc;
1081 const double one_m_Sr_loc = 1.0 - S_wr_loc;
1082 const double Se_trap_L1103 = 1.0 - S_gr.data()[mat_eN] / one_m_Sr_loc;
1085 const double z_cl = fmin(fmax(u_n, 1.0e-8), 1.0 - 1.0e-8);
1086 const double z_cl_old = fmin(fmax(u_n_old, 1.0e-8), 1.0 - 1.0e-8);
1087 const double p_cl = fmax(u_w_qp, 1.0e2);
1088 const double p_cl_old = fmax(u_w_qp_old, 1.0e2);
1094 const double N_old = fs_n_old.
rho_g*fs_n_old.
S_g
1095 + fs_n_old.
rho_a*(1.0 - fs_n_old.
S_g);
1096 const double m_n = phi_eN * N_cur * z_cl;
1097 const double m_n_old = phi_eN * N_old * z_cl_old;
1098 const double m_n_t = (m_n - m_n_old) / dt;
1102 const double S_a_qp = 1.0 - fs_n.
S_g;
1103 const double Se_a = fmin(fmax((S_a_qp - S_wr_loc) / one_m_Sr_loc, 0.0), 1.0);
1104 double KWr_a = 0.0, DKWr_a = 0.0, thW_a = 0.0, DthW_a = 0.0;
1105 double KNr_a = 0.0, DKNr_a = 0.0;
1106 double pc_a = 0.0, dpc_dSe_a = 0.0, d2pc_a = 0.0;
1109 thetaR.data()[mat_eN], thetaSR.data()[mat_eN], thW_a, DthW_a, KWr_a, DKWr_a);
1114 thetaR.data()[mat_eN], thetaSR.data()[mat_eN], thW_a, DthW_a, KWr_a, DKWr_a);
1118 KNr_a *= krn_end_eN;
1120 const double pcp_a = dpc_dSe_a / one_m_Sr_loc;
1126 const double rho_g_mass = fs_n.
rho_g*Mbar_g;
1127 const double rho_a_mass = fs_n.
rho_a*Mbar_a;
1133 for (
int I = 0; I < nSpace; I++) {
1134 double ua = 0.0, ug = 0.0;
1135 for (
int ii = a_rowptr.data()[I]; ii < a_rowptr.data()[I + 1]; ii++) {
1136 const int J = a_colind.data()[ii];
1137 const double gradSa_J = -(fs_n.
dS_g_dp*grad_u_w[J] + fs_n.
dS_g_dz*grad_u_n[J]);
1138 const double gp_a = grad_u_w[J] - rho_a_mass*gravity.data()[J];
1139 const double gp_g = grad_u_w[J] + pcp_a*gradSa_J - rho_g_mass*gravity.data()[J];
1140 ua -= (KWr_a*KWs_eN[ii]) * gp_a;
1141 ug -= (KNr_a*KWs_eN[ii]/mu_n) * gp_g;
1148 for (
int i = 0; i < nDOF_test_element; i++) {
1149 const double test_i = u_test_ref.data()[k * nDOF_test_element + i];
1151 elementResidual_n[i] += m_n_t * test_i * dV;
1152 lumped_w_n[i] += test_i * dV;
1153 for (
int I = 0; I < nSpace; I++) {
1154 elementResidual_n[i] -= F1[I] * u_grad_trial_qp[i * nSpace + I] * dV;
1159 for (
int i = 0; i < nDOF_test_element; i++) {
1160 const int gi = u_l2g.data()[eN * nDOF_test_element + i];
1161 elementResidual_n[i] -= injection_dof.data()[gi] * lumped_w_n[i];
1163 for (
int i = 0; i < nDOF_test_element; i++) {
1164 const int eN_i = eN * nDOF_test_element + i;
1165 globalResidual.data()[offset_n + stride_n * u_l2g.data()[eN_i]] += elementResidual_n[i];
1185 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++) {
1186 const int ebN = exteriorElementBoundariesArray.data()[ebNE];
1187 const int eN = elementBoundaryElementsArray.data()[ebN * 2 + 0];
1188 const int ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN * 2 + 0];
1189 const int eN_nDOF_trial_element = eN * nDOF_trial_element;
1190 const int mat_eN = elementMaterialTypes.data()[eN];
1191 const double phi_eN = thetaR.data()[mat_eN] + thetaSR.data()[mat_eN];
1192 const double alpha_eN = alpha.data()[mat_eN];
1193 const double krn_end_eN = krn_end.data()[mat_eN];
1194 const double n_vg_eN =
n.data()[mat_eN];
1195 const double *KWs_eN = &KWs.data()[mat_eN *
nnz];
1196 const double S_wr_loc = thetaR.data()[mat_eN] / phi_eN;
1197 const double one_m_Sr_loc = 1.0 - S_wr_loc;
1198 const double Se_trap_L1235 = 1.0 - S_gr.data()[mat_eN] / one_m_Sr_loc;
1200 double elementResidual_n_eb[nDOF_test_element];
1201 for (
int i = 0; i < nDOF_test_element; i++) elementResidual_n_eb[i] = 0.0;
1203 for (
int kb = 0; kb < nQuadraturePoints_elementBoundary; kb++) {
1204 const int ebNE_kb = ebNE * nQuadraturePoints_elementBoundary + kb;
1205 const int ebN_local_kb = ebN_local * nQuadraturePoints_elementBoundary + kb;
1206 const int ebN_local_kb_nSpace = ebN_local_kb * nSpace;
1208 double jac_ext[nSpace * nSpace], jacDet_ext, jacInv_ext[nSpace * nSpace];
1209 double boundaryJac_b[nSpace * (nSpace - 1)];
1210 double metricTensor_b[(nSpace - 1) * (nSpace - 1)];
1211 double metricTensorDetSqrt_b, dS_eb, normal_b[3];
1212 double xt_b, yt_b, zt_b, integralScaling_b;
1213 double x_eb, y_eb, z_eb;
1214 ck.calculateMapping_elementBoundary(eN, ebN_local, kb, ebN_local_kb,
1215 mesh_dof.data(), mesh_l2g.data(), mesh_trial_trace_ref.data(),
1216 mesh_grad_trial_trace_ref.data(), boundaryJac_ref.data(),
1217 jac_ext, jacDet_ext, jacInv_ext, boundaryJac_b, metricTensor_b,
1218 metricTensorDetSqrt_b, normal_ref.data(), normal_b,
1220 ck.calculateMappingVelocity_elementBoundary(eN, ebN_local, kb, ebN_local_kb,
1221 mesh_velocity_dof.data(), mesh_l2g.data(), mesh_trial_trace_ref.data(),
1222 xt_b, yt_b, zt_b, normal_b, boundaryJac_b, metricTensor_b,
1224 dS_eb = ((1.0 - MOVING_DOMAIN) * metricTensorDetSqrt_b
1225 + MOVING_DOMAIN * integralScaling_b) * dS_ref.data()[kb];
1227 double u_grad_trial_trace_b[nDOF_trial_element * nSpace];
1228 ck.gradTrialFromRef(
1229 &u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace * nDOF_trial_element],
1230 jacInv_ext, u_grad_trial_trace_b);
1231 double u_w_ext_b = 0.0, u_n_ext_b = 0.0;
1232 double grad_u_w_ext_b[nSpace], grad_u_n_ext_b[nSpace];
1233 ck.valFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element],
1234 &u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element],
1236 ck.valFromDOF(u_dof_n.data(), &u_l2g.data()[eN_nDOF_trial_element],
1237 &u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element],
1239 ck.gradFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element],
1240 u_grad_trial_trace_b, grad_u_w_ext_b);
1241 ck.gradFromDOF(u_dof_n.data(), &u_l2g.data()[eN_nDOF_trial_element],
1242 u_grad_trial_trace_b, grad_u_n_ext_b);
1244 const int isDir_n = isDOFBoundary_n.data()[ebNE_kb];
1245 const double bc_u_n_ext_b = isDir_n * ebqe_bc_u_n_ext.data()[ebNE_kb]
1246 + (1 - isDir_n) * u_n_ext_b;
1259 const double z_clb = fmin(fmax(u_n_ext_b, 1.0e-8), 1.0 - 1.0e-8);
1260 const double p_clb = fmax(u_w_ext_b, 1.0e2);
1263 const double Se_ab = fmin(fmax((1.0 - fsb.
S_g - S_wr_loc)/one_m_Sr_loc, 0.0), 1.0);
1264 double KWr_b=0,DKWr_b=0,thW_b=0,DthW_b=0,KNr_b=0,DKNr_b=0,pc_b=0,dpc_dSe_b=0,d2pc_b=0;
1274 KNr_b *= krn_end_eN;
1275 const double pcpb = dpc_dSe_b / one_m_Sr_loc;
1278 const double rho_g_mass_b = fsb.
rho_g*Mbar_gb;
1279 const double rho_a_mass_b = fsb.
rho_a*Mbar_ab;
1280 double F_n_dot_n = 0.0;
1281 for (
int I = 0; I < nSpace; I++) {
1282 double ua = 0.0, ug = 0.0;
1283 for (
int ii = a_rowptr.data()[I]; ii < a_rowptr.data()[I + 1]; ii++) {
1284 const int J = a_colind.data()[ii];
1285 const double gradSa_J = -(fsb.
dS_g_dp*grad_u_w_ext_b[J] + fsb.
dS_g_dz*grad_u_n_ext_b[J]);
1286 const double gp_a = grad_u_w_ext_b[J] - rho_a_mass_b*gravity.data()[J];
1287 const double gp_g = grad_u_w_ext_b[J] + pcpb*gradSa_J - rho_g_mass_b*gravity.data()[J];
1288 ua -= (KWr_b*KWs_eN[ii]) * gp_a;
1289 ug -= (KNr_b*KWs_eN[ii]/mu_n) * gp_g;
1291 F_n_dot_n += (fsb.
rho_g*fsb.
Y*ug + fsb.
rho_a*fsb.
X*ua) * normal_b[I];
1300 double Kw_rep = 0.0;
1301 for (
int ii = 0; ii <
nnz; ii++) Kw_rep = fmax(Kw_rep, fabs(KWs_eN[ii]));
1302 const double a_n_scale = (fsb.
rho_g*fsb.
Y*KNr_b/mu_n
1303 + fsb.
rho_a*fsb.
X*KWr_b) * Kw_rep;
1304 const double pen_eff = ebqe_penalty_ext.data()[ebNE_kb] * a_n_scale;
1305 F_n_dot_n += pen_eff * (u_n_ext_b - bc_u_n_ext_b);
1310 for (
int i = 0; i < nDOF_test_element; i++) {
1311 const double test_i_dS = u_test_trace_ref.data()[
1312 ebN_local_kb * nDOF_test_element + i] * dS_eb;
1313 elementResidual_n_eb[i] += F_n_dot_n * test_i_dS;
1317 for (
int i = 0; i < nDOF_test_element; i++) {
1318 const int gi = u_l2g.data()[eN * nDOF_test_element + i];
1319 globalResidual.data()[offset_n + stride_n * gi] += elementResidual_n_eb[i];
1326 xt::pyarray<double> &mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
1327 xt::pyarray<double> &mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
1328 xt::pyarray<double> &mesh_dof = args.
array<
double>(
"mesh_dof");
1329 xt::pyarray<double> &mesh_velocity_dof = args.
array<
double>(
"mesh_velocity_dof");
1330 double MOVING_DOMAIN = args.
scalar<
double>(
"MOVING_DOMAIN");
1331 xt::pyarray<int> &mesh_l2g = args.
array<
int>(
"mesh_l2g");
1332 xt::pyarray<double> &dV_ref = args.
array<
double>(
"dV_ref");
1333 xt::pyarray<double> &u_trial_ref = args.
array<
double>(
"u_trial_ref");
1334 xt::pyarray<double> &u_grad_trial_ref = args.
array<
double>(
"u_grad_trial_ref");
1335 xt::pyarray<double> &u_test_ref = args.
array<
double>(
"u_test_ref");
1336 xt::pyarray<double> &u_grad_test_ref = args.
array<
double>(
"u_grad_test_ref");
1337 xt::pyarray<double> &mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
1338 xt::pyarray<double> &mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
1339 xt::pyarray<double> &dS_ref = args.
array<
double>(
"dS_ref");
1340 xt::pyarray<double> &u_trial_trace_ref = args.
array<
double>(
"u_trial_trace_ref");
1341 xt::pyarray<double> &u_grad_trial_trace_ref = args.
array<
double>(
"u_grad_trial_trace_ref");
1342 xt::pyarray<double> &u_test_trace_ref = args.
array<
double>(
"u_test_trace_ref");
1343 xt::pyarray<double> &u_grad_test_trace_ref = args.
array<
double>(
"u_grad_test_trace_ref");
1344 xt::pyarray<double> &normal_ref = args.
array<
double>(
"normal_ref");
1345 xt::pyarray<double> &boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
1346 int nElements_global = args.
scalar<
int>(
"nElements_global");
1347 xt::pyarray<double> &ebqe_penalty_ext = args.
array<
double>(
"ebqe_penalty_ext");
1348 xt::pyarray<int> &elementMaterialTypes = args.
array<
int>(
"elementMaterialTypes");
1349 xt::pyarray<int> &isSeepageFace = args.
array<
int>(
"isSeepageFace");
1350 xt::pyarray<int> &a_rowptr = args.
array<
int>(
"a_rowptr");
1351 xt::pyarray<int> &a_colind = args.
array<
int>(
"a_colind");
1352 double rho = args.
scalar<
double>(
"rho");
1353 double beta = args.
scalar<
double>(
"beta");
1356 xt::pyarray<double> &q_rho = args.
array<
double>(
"q_rho");
1357 xt::pyarray<double> &ebqe_rho = args.
array<
double>(
"ebqe_rho");
1365 xt::pyarray<double> &c_dof_jac = args.
array<
double>(
"c_dof");
1366 const double k_d_jac = args.
scalar<
double>(
"k_d");
1367 const double c_sat_jac = args.
scalar<
double>(
"c_sat");
1369 xt::pyarray<double> &gravity = args.
array<
double>(
"gravity");
1370 xt::pyarray<double> &alpha = args.
array<
double>(
"alpha");
1371 xt::pyarray<double> &
n = args.
array<
double>(
"n");
1372 xt::pyarray<double> &thetaR = args.
array<
double>(
"thetaR");
1373 xt::pyarray<double> &thetaSR = args.
array<
double>(
"thetaSR");
1374 xt::pyarray<double> &KWs = args.
array<
double>(
"KWs");
1375 xt::pyarray<double> &krn_end = args.
array<
double>(
"krn_end");
1376 xt::pyarray<double> &S_gr = args.
array<
double>(
"S_gr");
1377 double mu_n = args.
scalar<
double>(
"mu_n");
1378 double useMetrics = args.
scalar<
double>(
"useMetrics");
1379 double alphaBDF = args.
scalar<
double>(
"alphaBDF");
1380 int lag_shockCapturing = args.
scalar<
int>(
"lag_shockCapturing");
1381 double shockCapturingDiffusion = args.
scalar<
double>(
"shockCapturingDiffusion");
1383 double VMS = args.
scalar<
double>(
"VMS");
1384 xt::pyarray<int> &u_l2g = args.
array<
int>(
"u_l2g");
1385 xt::pyarray<double> &elementDiameter = args.
array<
double>(
"elementDiameter");
1386 xt::pyarray<double> &u_dof = args.
array<
double>(
"u_dof");
1389 xt::pyarray<double> &u_dof_n = args.
array<
double>(
"u_dof_n");
1390 xt::pyarray<double> &velocity = args.
array<
double>(
"velocity");
1391 xt::pyarray<double> &q_m_betaBDF = args.
array<
double>(
"q_m_betaBDF");
1392 xt::pyarray<double> &cfl = args.
array<
double>(
"cfl");
1393 xt::pyarray<double> &q_numDiff_u = args.
array<
double>(
"q_numDiff_u");
1394 xt::pyarray<double> &q_numDiff_u_last = args.
array<
double>(
"q_numDiff_u_last");
1395 xt::pyarray<int> &csrRowIndeces_u_u = args.
array<
int>(
"csrRowIndeces_u_u");
1396 xt::pyarray<int> &csrColumnOffsets_u_u = args.
array<
int>(
"csrColumnOffsets_u_u");
1397 xt::pyarray<double> &globalJacobian = args.
array<
double>(
"globalJacobian");
1402 const double dt_n = args.
scalar<
double>(
"dt");
1406 const double rho_n = args.
scalar<
double>(
"rho_n");
1407 const double p_ref_n = args.
scalar<
double>(
"p_ref_n");
1408 const bool rho_n_compressible = (p_ref_n > 0.0);
1409 const double c_n = rho_n_compressible ? (rho_n / p_ref_n) : 0.0;
1410 xt::pyarray<int> &csrRowIndeces_n_n = args.
array<
int>(
"csrRowIndeces_n_n");
1414 xt::pyarray<int> &csrRowIndeces_n_w = args.
array<
int>(
"csrRowIndeces_n_w");
1415 xt::pyarray<int> &csrColumnOffsets_n_n = args.
array<
int>(
"csrColumnOffsets_n_n");
1416 xt::pyarray<int> &csrColumnOffsets_n_w = args.
array<
int>(
"csrColumnOffsets_n_w");
1420 xt::pyarray<int> &csrColumnOffsets_eb_n_n = args.
array<
int>(
"csrColumnOffsets_eb_n_n");
1421 xt::pyarray<int> &csrColumnOffsets_eb_n_w = args.
array<
int>(
"csrColumnOffsets_eb_n_w");
1427 xt::pyarray<int> &csrRowIndeces_w_n = args.
array<
int>(
"csrRowIndeces_w_n");
1428 xt::pyarray<int> &csrColumnOffsets_w_n = args.
array<
int>(
"csrColumnOffsets_w_n");
1431 xt::pyarray<int> &csrColumnOffsets_eb_w_n = args.
array<
int>(
"csrColumnOffsets_eb_w_n");
1436 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
1437 xt::pyarray<int> &exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
1438 xt::pyarray<int> &elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
1439 xt::pyarray<int> &elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
1440 xt::pyarray<double> &ebqe_velocity_ext = args.
array<
double>(
"ebqe_velocity_ext");
1441 xt::pyarray<int> &isDOFBoundary_u = args.
array<
int>(
"isDOFBoundary_u");
1442 xt::pyarray<double> &ebqe_bc_u_ext = args.
array<
double>(
"ebqe_bc_u_ext");
1444 xt::pyarray<int> &isDOFBoundary_n = args.
array<
int>(
"isDOFBoundary_n");
1445 xt::pyarray<double> &ebqe_bc_u_n_ext = args.
array<
double>(
"ebqe_bc_u_n_ext");
1446 xt::pyarray<int> &isFluxBoundary_u = args.
array<
int>(
"isFluxBoundary_u");
1447 xt::pyarray<double> &ebqe_bc_flux_ext = args.
array<
double>(
"ebqe_bc_flux_ext");
1448 xt::pyarray<int> &csrColumnOffsets_eb_u_u = args.
array<
int>(
"csrColumnOffsets_eb_u_u");
1449 int LUMPED_MASS_MATRIX = args.
scalar<
int>(
"LUMPED_MASS_MATRIX");
1450 assert(a_rowptr.data()[nSpace] ==
nnz);
1451 assert(a_rowptr.data()[nSpace] == nSpace);
1452 double Ct_sge = 4.0;
1457 for (
int eN = 0; eN < nElements_global; eN++) {
1458 double elementJacobian_u_u[nDOF_test_element][nDOF_trial_element];
1460 double elementJacobian_u_n[nDOF_test_element][nDOF_trial_element];
1461 for (
int i = 0; i < nDOF_test_element; i++) {
1462 for (
int j = 0; j < nDOF_trial_element; j++) {
1463 elementJacobian_u_u[i][j] = 0.0;
1464 elementJacobian_u_n[i][j] = 0.0;
1467 for (
int k = 0; k < nQuadraturePoints_element; k++) {
1468 int eN_k = eN * nQuadraturePoints_element + k,
1469 eN_k_nSpace = eN_k * nSpace,
1470 eN_nDOF_trial_element = eN * nDOF_trial_element;
1473 double u = 0.0, grad_u[nSpace], m = 0.0, dm = 0.0,
f[nSpace],
df[nSpace], a[
nnz], da[
nnz], as[
nnz], m_t = 0.0, dm_t = 0.0, dpdeResidual_u_u[nDOF_trial_element], Lstar_u[nDOF_test_element], dsubgridError_u_u[nDOF_trial_element], tau = 0.0, tau0 = 0.0, tau1 = 0.0, jac[nSpace * nSpace], jacDet, jacInv[nSpace * nSpace], u_grad_trial[nDOF_trial_element * nSpace], dV, u_test_dV[nDOF_test_element], u_grad_test_dV[nDOF_test_element * nSpace], x, y,
z,
xt, yt, zt, G[nSpace * nSpace], G_dd_G, tr_G;
1478 ck.calculateMapping_element(eN, k, mesh_dof.data(), mesh_l2g.data(), mesh_trial_ref.data(), mesh_grad_trial_ref.data(), jac, jacDet, jacInv, x, y,
z);
1479 ck.calculateMappingVelocity_element(eN, k, mesh_velocity_dof.data(), mesh_l2g.data(), mesh_trial_ref.data(),
xt, yt, zt);
1481 dV = fabs(jacDet) * dV_ref.data()[k];
1482 ck.calculateG(jacInv, G, G_dd_G, tr_G);
1484 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k * nDOF_trial_element * nSpace], jacInv, u_grad_trial);
1486 ck.valFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], &u_trial_ref.data()[k * nDOF_trial_element],
u);
1488 ck.gradFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], u_grad_trial, grad_u);
1490 for (
int j = 0; j < nDOF_trial_element; j++) {
1491 u_test_dV[j] = u_test_ref.data()[k * nDOF_trial_element + j] * dV;
1492 for (
int I = 0; I < nSpace; I++) {
1493 u_grad_test_dV[j * nSpace + I] = u_grad_trial[j * nSpace + I] * dV;
1504 const int mat_eN0 = elementMaterialTypes.data()[eN];
1505 const double alpha_eN0 = alpha.data()[mat_eN0];
1506 const double n_vg_eN0 =
n.data()[mat_eN0];
1507 const double krn_end_eN0 = krn_end.data()[mat_eN0];
1508 const double *KWs_eN0 = &KWs.data()[mat_eN0 *
nnz];
1509 const double phi_eN0 = thetaR.data()[mat_eN0] + thetaSR.data()[mat_eN0];
1510 const double S_wr_loc0 = thetaR.data()[mat_eN0] / phi_eN0;
1511 const double one_m_Sr0 = 1.0 - S_wr_loc0;
1512 const double Se_trap_L1547 = 1.0 - S_gr.data()[mat_eN0] / one_m_Sr0;
1513 double u_n_qp = 0.0, grad_u_n[nSpace];
1514 ck.valFromDOF(u_dof_n.data(), &u_l2g.data()[eN_nDOF_trial_element],
1515 &u_trial_ref.data()[k * nDOF_trial_element], u_n_qp);
1516 ck.gradFromDOF(u_dof_n.data(), &u_l2g.data()[eN_nDOF_trial_element],
1517 u_grad_trial, grad_u_n);
1518 const double z_cl0 = fmin(fmax(u_n_qp, 1.0e-8), 1.0 - 1.0e-8);
1519 const double p_cl0 = fmax(
u, 1.0e2);
1522 const double S_g0 = fs0.
S_g, Sa0 = 1.0 - S_g0;
1524 const double N0 = fs0.
rho_g*S_g0 + fs0.
rho_a*Sa0;
1529 const double dm0_dp = phi_eN0 * dN0_dp * (1.0 - z_cl0);
1530 const double dm0_dz = phi_eN0 * (dN0_dz * (1.0 - z_cl0) - N0);
1532 const double Se_raw0 = (Sa0 - S_wr_loc0) / one_m_Sr0;
1533 double Se_a0, dSe0_dp, dSe0_dz;
1534 if (Se_raw0 <= 0.0) { Se_a0 = 0.0; dSe0_dp = 0.0; dSe0_dz = 0.0; }
1535 else if (Se_raw0 >= 1.0) { Se_a0 = 1.0; dSe0_dp = 0.0; dSe0_dz = 0.0; }
1536 else { Se_a0 = Se_raw0; dSe0_dp = -fs0.
dS_g_dp/one_m_Sr0; dSe0_dz = -fs0.
dS_g_dz/one_m_Sr0; }
1537 double KWr0=0.0, DKWr0=0.0, thW0=0.0, DthW0=0.0, KNr0=0.0, DKNr0=0.0;
1538 double pc0=0.0, dpc_dSe0=0.0, d2pc0=0.0;
1541 thetaR.data()[mat_eN0], thetaSR.data()[mat_eN0], thW0, DthW0, KWr0, DKWr0);
1546 thetaR.data()[mat_eN0], thetaSR.data()[mat_eN0], thW0, DthW0, KWr0, DKWr0);
1550 KNr0 *= krn_end_eN0; DKNr0 *= krn_end_eN0;
1551 const double pcp0 = dpc_dSe0 / one_m_Sr0;
1552 const double dpcp0_dp = (d2pc0 / one_m_Sr0) * dSe0_dp;
1553 const double dpcp0_dz = (d2pc0 / one_m_Sr0) * dSe0_dz;
1557 const double rho_g_mass0 = fs0.
rho_g*Mbar_g0, rho_a_mass0 = fs0.
rho_a*Mbar_a0;
1559 const double drgm0_dz = fs0.
rho_g*fs0.
dY_dz*dMm0;
1563 const double Ag0 = fs0.
rho_g*(1.0-fs0.
Y), Aa0 = fs0.
rho_a*(1.0-fs0.
X);
1565 const double dAg0_dz = - fs0.
rho_g*fs0.
dY_dz;
1569 double ug0[nSpace], ua0[nSpace];
1570 double dug0_dp[nSpace], dug0_dz[nSpace], dua0_dp[nSpace], dua0_dz[nSpace];
1571 for (
int I = 0; I < nSpace; I++) {
1572 double ugI=0.0, uaI=0.0, dugp=0.0, dugz=0.0, duap=0.0, duaz=0.0;
1573 for (
int ii = a_rowptr.data()[I]; ii < a_rowptr.data()[I + 1]; ii++) {
1574 const int J = a_colind.data()[ii];
1575 const double Kii = KWs_eN0[ii];
1576 const double Mob_g = KNr0*Kii/mu_n, Mob_a = KWr0*Kii;
1577 const double dMobg_dp = (DKNr0*Kii/mu_n)*dSe0_dp, dMobg_dz = (DKNr0*Kii/mu_n)*dSe0_dz;
1578 const double dMoba_dp = (DKWr0*Kii)*dSe0_dp, dMoba_dz = (DKWr0*Kii)*dSe0_dz;
1579 const double gJ = gravity.data()[J];
1580 const double gradSa = -(fs0.
dS_g_dp*grad_u[J] + fs0.
dS_g_dz*grad_u_n[J]);
1581 const double gp_a = grad_u[J] - rho_a_mass0*gJ;
1582 const double gp_g = grad_u[J] + pcp0*gradSa - rho_g_mass0*gJ;
1583 ugI -= Mob_g*gp_g; uaI -= Mob_a*gp_a;
1586 const double dgpg_dp = dpcp0_dp*gradSa + pcp0*dgradSa_dp - drgm0_dp*gJ;
1587 const double dgpg_dz = dpcp0_dz*gradSa + pcp0*dgradSa_dz - drgm0_dz*gJ;
1588 dugp -= dMobg_dp*gp_g + Mob_g*dgpg_dp;
1589 dugz -= dMobg_dz*gp_g + Mob_g*dgpg_dz;
1590 duap -= dMoba_dp*gp_a + Mob_a*(-dram0_dp*gJ);
1591 duaz -= dMoba_dz*gp_a + Mob_a*(-dram0_dz*gJ);
1593 ug0[I]=ugI; ua0[I]=uaI;
1594 dug0_dp[I]=dugp; dug0_dz[I]=dugz; dua0_dp[I]=duap; dua0_dz[I]=duaz;
1600 for (
int i = 0; i < nDOF_test_element; i++) {
1601 const int i_nSpace = i * nSpace;
1602 double Sval_p = 0.0, Sval_z = 0.0;
1603 for (
int I = 0; I < nSpace; I++) {
1604 const double gNiI = u_grad_trial[i_nSpace + I];
1605 Sval_p += (dAg0_dp*ug0[I] + Ag0*dug0_dp[I] + dAa0_dp*ua0[I] + Aa0*dua0_dp[I]) * gNiI;
1606 Sval_z += (dAg0_dz*ug0[I] + Ag0*dug0_dz[I] + dAa0_dz*ua0[I] + Aa0*dua0_dz[I]) * gNiI;
1608 for (
int j = 0; j < nDOF_trial_element; j++) {
1609 const int j_nSpace = j * nSpace;
1610 const double trial_j = u_trial_ref.data()[k * nDOF_trial_element + j];
1611 double Sgrad_p = 0.0, Sgrad_z = 0.0;
1612 for (
int I = 0; I < nSpace; I++) {
1613 const double gNiI = u_grad_trial[i_nSpace + I];
1614 for (
int ii = a_rowptr.data()[I]; ii < a_rowptr.data()[I + 1]; ii++) {
1615 const int J = a_colind.data()[ii];
1616 const double Kii = KWs_eN0[ii];
1617 const double Mob_g = KNr0*Kii/mu_n, Mob_a = KWr0*Kii;
1618 const double gNjJ = u_grad_trial[j_nSpace + J];
1619 const double dFdgp = Ag0*(-Mob_g*(1.0 - pcp0*fs0.
dS_g_dp)) + Aa0*(-Mob_a);
1620 const double dFdgz = Ag0*(-Mob_g*(-pcp0*fs0.
dS_g_dz));
1621 Sgrad_p += dFdgp * gNjJ * gNiI;
1622 Sgrad_z += dFdgz * gNjJ * gNiI;
1625 elementJacobian_u_u[i][j] += alphaBDF * dm0_dp * trial_j * u_test_dV[i];
1626 elementJacobian_u_u[i][j] -= (Sval_p * trial_j + Sgrad_p) * dV;
1627 elementJacobian_u_u[i][j] +=
VMS *
ck.NumericalDiffusionJacobian(
1628 q_numDiff_u_last[eN_k], &u_grad_trial[j_nSpace], &u_grad_test_dV[i_nSpace]);
1629 elementJacobian_u_n[i][j] += alphaBDF * dm0_dz * trial_j * u_test_dV[i];
1630 elementJacobian_u_n[i][j] -= (Sval_z * trial_j + Sgrad_z) * dV;
1637 for (
int i = 0; i < nDOF_test_element; i++) {
1638 int eN_i = eN * nDOF_test_element + i;
1639 for (
int j = 0; j < nDOF_trial_element; j++) {
1640 int eN_i_j = eN_i * nDOF_trial_element + j;
1641 globalJacobian.data()[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_u_u[eN_i_j]] += elementJacobian_u_u[i][j];
1644 globalJacobian.data()[csrRowIndeces_w_n.data()[eN_i] + csrColumnOffsets_w_n.data()[eN_i_j]] += elementJacobian_u_n[i][j];
1651 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++) {
1652 int ebN = exteriorElementBoundariesArray.data()[ebNE];
1653 int eN = elementBoundaryElementsArray.data()[ebN * 2 + 0], ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN * 2 + 0], eN_nDOF_trial_element = eN * nDOF_trial_element;
1654 for (
int kb = 0; kb < nQuadraturePoints_elementBoundary; kb++) {
1655 int ebNE_kb = ebNE * nQuadraturePoints_elementBoundary + kb, ebNE_kb_nSpace = ebNE_kb * nSpace, ebN_local_kb = ebN_local * nQuadraturePoints_elementBoundary + kb, ebN_local_kb_nSpace = ebN_local_kb * nSpace;
1657 double u_ext = 0.0, grad_u_ext[nSpace], m_ext = 0.0, dm_ext = 0.0, f_ext[nSpace], df_ext[nSpace], a_ext[
nnz], da_ext[
nnz], as_ext[
nnz], dflux_u_u_ext = 0.0, bc_u_ext = 0.0,
1659 bc_m_ext = 0.0, bc_dm_ext = 0.0, bc_f_ext[nSpace], bc_df_ext[nSpace], bc_a_ext[
nnz], bc_da_ext[
nnz], bc_as_ext[
nnz], fluxJacobian_u_u[nDOF_trial_element], jac_ext[nSpace * nSpace], jacDet_ext, jacInv_ext[nSpace * nSpace], boundaryJac[nSpace * (nSpace - 1)], metricTensor[(nSpace - 1) * (nSpace - 1)], metricTensorDetSqrt, dS, u_test_dS[nDOF_test_element], u_grad_trial_trace[nDOF_trial_element * nSpace], normal[3], x_ext, y_ext, z_ext, xt_ext, yt_ext, zt_ext, integralScaling, G[nSpace * nSpace], G_dd_G, tr_G;
1663 ck.calculateMapping_elementBoundary(eN, ebN_local, kb, ebN_local_kb, mesh_dof.data(), mesh_l2g.data(), mesh_trial_trace_ref.data(), mesh_grad_trial_trace_ref.data(), boundaryJac_ref.data(), jac_ext, jacDet_ext, jacInv_ext, boundaryJac, metricTensor, metricTensorDetSqrt,
1664 normal_ref.data(), normal, x_ext, y_ext, z_ext);
1665 ck.calculateMappingVelocity_elementBoundary(eN, ebN_local, kb, ebN_local_kb, mesh_velocity_dof.data(), mesh_l2g.data(), mesh_trial_trace_ref.data(), xt_ext, yt_ext, zt_ext, normal, boundaryJac, metricTensor, integralScaling);
1666 dS = ((1.0 - MOVING_DOMAIN) * metricTensorDetSqrt + MOVING_DOMAIN * integralScaling) * dS_ref.data()[kb];
1667 ck.calculateG(jacInv_ext, G, G_dd_G, tr_G);
1670 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace * nDOF_trial_element], jacInv_ext, u_grad_trial_trace);
1672 ck.valFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], &u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element], u_ext);
1673 ck.gradFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], u_grad_trial_trace, grad_u_ext);
1675 for (
int j = 0; j < nDOF_trial_element; j++) { u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb * nDOF_test_element + j] * dS; }
1679 bc_u_ext = isDOFBoundary_u.data()[ebNE_kb] * ebqe_bc_u_ext.data()[ebNE_kb] + (1 - isDOFBoundary_u.data()[ebNE_kb]) * u_ext;
1683 double Kr, dKr, thetaW, thetaW_bc;
1684 const double rho_ext = ebqe_rho.data()[ebNE_kb];
1687 double dm_du_n_ext = 0.0, dkr_du_n_ext = 0.0;
1688 double df_du_n_ext[nSpace];
1689 double da_du_n_ext[
nnz];
1690 double bc_dm_du_n = 0.0, bc_dkr_du_n = 0.0;
1691 double bc_df_du_n[nSpace];
1692 double bc_da_du_n[
nnz];
1693 for (
int I = 0; I < nSpace; I++) { df_du_n_ext[I] = 0.0; bc_df_du_n[I] = 0.0; }
1694 for (
int ii = 0; ii <
nnz; ii++) { da_du_n_ext[ii] = 0.0; bc_da_du_n[ii] = 0.0; }
1697 double u_n_ext_qp_outer = 0.0;
1698 double bc_u_n_ext_qp_outer = 0.0;
1700 double u_n_ext_qp = 0.0;
1701 ck.valFromDOF(u_dof_n.data(), &u_l2g.data()[eN_nDOF_trial_element],
1702 &u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element], u_n_ext_qp);
1703 const double bc_u_n_ext_qp = isDOFBoundary_n.data()[ebNE_kb] * ebqe_bc_u_n_ext.data()[ebNE_kb]
1704 + (1 - isDOFBoundary_n.data()[ebNE_kb]) * u_n_ext_qp;
1705 u_n_ext_qp_outer = u_n_ext_qp;
1706 bc_u_n_ext_qp_outer = bc_u_n_ext_qp;
1708 alpha.data()[elementMaterialTypes.data()[eN]],
n.data()[elementMaterialTypes.data()[eN]],
1709 thetaR.data()[elementMaterialTypes.data()[eN]], thetaSR.data()[elementMaterialTypes.data()[eN]],
1710 &KWs.data()[elementMaterialTypes.data()[eN] *
nnz], u_ext, u_n_ext_qp,
1711 m_ext, dm_ext, dm_du_n_ext, f_ext, df_ext, df_du_n_ext, a_ext, da_ext, da_du_n_ext,
1712 as_ext, Kr, dKr, dkr_du_n_ext, thetaW);
1714 alpha.data()[elementMaterialTypes.data()[eN]],
n.data()[elementMaterialTypes.data()[eN]],
1715 thetaR.data()[elementMaterialTypes.data()[eN]], thetaSR.data()[elementMaterialTypes.data()[eN]],
1716 &KWs.data()[elementMaterialTypes.data()[eN] *
nnz], bc_u_ext, bc_u_n_ext_qp,
1717 bc_m_ext, bc_dm_ext, bc_dm_du_n, bc_f_ext, bc_df_ext, bc_df_du_n, bc_a_ext, bc_da_ext, bc_da_du_n,
1718 bc_as_ext, Kr, dKr, bc_dkr_du_n, thetaW_bc);
1727 double fluxJacobian_u_n[nDOF_trial_element];
1728 for (
int j = 0; j < nDOF_trial_element; j++) { fluxJacobian_u_u[j] = 0.0; fluxJacobian_u_n[j] = 0.0; }
1729 if (isDOFBoundary_u.data()[ebNE_kb]) {
1730 double grad_u_n_ext[nSpace];
1731 ck.gradFromDOF(u_dof_n.data(), &u_l2g.data()[eN_nDOF_trial_element],
1732 u_grad_trial_trace, grad_u_n_ext);
1733 const int mat_b = elementMaterialTypes.data()[eN];
1734 const double alpha_b = alpha.data()[mat_b];
1735 const double n_vg_b =
n.data()[mat_b];
1736 const double krn_end_b = krn_end.data()[mat_b];
1737 const double *KWs_b = &KWs.data()[mat_b *
nnz];
1738 const double phi_b = thetaR.data()[mat_b] + thetaSR.data()[mat_b];
1739 const double S_wr_b = thetaR.data()[mat_b] / phi_b;
1740 const double one_m_Sr_b = 1.0 - S_wr_b;
1741 const double Se_trap_L1775 = 1.0 - S_gr.data()[mat_b] / one_m_Sr_b;
1742 const double z_clb = fmin(fmax(u_n_ext_qp_outer, 1.0e-8), 1.0 - 1.0e-8);
1743 const double p_clb = fmax(u_ext, 1.0e2);
1746 const double Sab = 1.0 - fsb.
S_g;
1747 const double Se_rawb = (Sab - S_wr_b)/one_m_Sr_b;
1748 double Se_ab, dSeb_dp, dSeb_dz;
1749 if (Se_rawb<=0.0){Se_ab=0.0;dSeb_dp=0.0;dSeb_dz=0.0;}
1750 else if (Se_rawb>=1.0){Se_ab=1.0;dSeb_dp=0.0;dSeb_dz=0.0;}
1751 else {Se_ab=Se_rawb;dSeb_dp=-fsb.
dS_g_dp/one_m_Sr_b;dSeb_dz=-fsb.
dS_g_dz/one_m_Sr_b;}
1752 double KWrb=0,DKWrb=0,thWb=0,DthWb=0,KNrb=0,DKNrb=0,pcb=0,dpc_dSeb=0,d2pcb=0;
1762 KNrb *= krn_end_b; DKNrb *= krn_end_b;
1763 const double pcpb = dpc_dSeb / one_m_Sr_b;
1764 const double dpcpb_dp = (d2pcb / one_m_Sr_b) * dSeb_dp;
1765 const double dpcpb_dz = (d2pcb / one_m_Sr_b) * dSeb_dz;
1769 const double rgmb = fsb.
rho_g*Mbar_gb, ramb = fsb.
rho_a*Mbar_ab;
1771 const double drgmb_dz = fsb.
rho_g*fsb.
dY_dz*dMmb;
1774 const double Agb = fsb.
rho_g*(1.0-fsb.
Y), Aab = fsb.
rho_a*(1.0-fsb.
X);
1776 const double dAgb_dz = - fsb.
rho_g*fsb.
dY_dz;
1779 double ugb[nSpace], uab[nSpace], dugb_dp[nSpace], dugb_dz[nSpace], duab_dp[nSpace], duab_dz[nSpace];
1780 for (
int I=0;I<nSpace;I++){
1781 double ugI=0,uaI=0,dugp=0,dugz=0,duap=0,duaz=0;
1782 for (
int ii=a_rowptr.data()[I];ii<a_rowptr.data()[I+1];ii++){
1783 const int J=a_colind.data()[ii];
1784 const double Kii=KWs_b[ii];
1785 const double Mob_g=KNrb*Kii/mu_n, Mob_a=KWrb*Kii;
1786 const double dMobg_dp=(DKNrb*Kii/mu_n)*dSeb_dp, dMobg_dz=(DKNrb*Kii/mu_n)*dSeb_dz;
1787 const double dMoba_dp=(DKWrb*Kii)*dSeb_dp, dMoba_dz=(DKWrb*Kii)*dSeb_dz;
1788 const double gJ=gravity.data()[J];
1789 const double gradSa=-(fsb.
dS_g_dp*grad_u_ext[J]+fsb.
dS_g_dz*grad_u_n_ext[J]);
1790 const double gp_a=grad_u_ext[J]-ramb*gJ;
1791 const double gp_g=grad_u_ext[J]+pcpb*gradSa-rgmb*gJ;
1792 ugI-=Mob_g*gp_g; uaI-=Mob_a*gp_a;
1795 const double dgpg_dp=dpcpb_dp*gradSa+pcpb*dgradSa_dp-drgmb_dp*gJ;
1796 const double dgpg_dz=dpcpb_dz*gradSa+pcpb*dgradSa_dz-drgmb_dz*gJ;
1797 dugp-=dMobg_dp*gp_g+Mob_g*dgpg_dp;
1798 dugz-=dMobg_dz*gp_g+Mob_g*dgpg_dz;
1799 duap-=dMoba_dp*gp_a+Mob_a*(-dramb_dp*gJ);
1800 duaz-=dMoba_dz*gp_a+Mob_a*(-dramb_dz*gJ);
1802 ugb[I]=ugI;uab[I]=uaI;dugb_dp[I]=dugp;dugb_dz[I]=dugz;duab_dp[I]=duap;duab_dz[I]=duaz;
1804 const double penb = ebqe_penalty_ext.data()[ebNE_kb];
1805 for (
int j=0;j<nDOF_trial_element;j++){
1806 const int j_nSpace = j*nSpace;
1807 const double trial_j = u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element + j];
1808 double accp=0.0, accz=0.0;
1809 for (
int I=0;I<nSpace;I++){
1810 const double nI = normal[I];
1811 const double valp = dAgb_dp*ugb[I]+Agb*dugb_dp[I]+dAab_dp*uab[I]+Aab*duab_dp[I];
1812 const double valz = dAgb_dz*ugb[I]+Agb*dugb_dz[I]+dAab_dz*uab[I]+Aab*duab_dz[I];
1813 double gradp=0.0, gradz=0.0;
1814 for (
int ii=a_rowptr.data()[I];ii<a_rowptr.data()[I+1];ii++){
1815 const int J=a_colind.data()[ii];
1816 const double Kii=KWs_b[ii];
1817 const double Mob_g=KNrb*Kii/mu_n, Mob_a=KWrb*Kii;
1818 const double gNjJ = u_grad_trial_trace[j_nSpace + J];
1819 const double dFdgp = Agb*(-Mob_g*(1.0 - pcpb*fsb.
dS_g_dp)) + Aab*(-Mob_a);
1820 const double dFdgz = Agb*(-Mob_g*(-pcpb*fsb.
dS_g_dz));
1821 gradp += dFdgp*gNjJ;
1822 gradz += dFdgz*gNjJ;
1824 accp += nI*(valp*trial_j + gradp);
1825 accz += nI*(valz*trial_j + gradz);
1827 fluxJacobian_u_u[j] = accp + penb*trial_j;
1828 fluxJacobian_u_n[j] = accz;
1834 for (
int i = 0; i < nDOF_test_element; i++) {
1835 int eN_i = eN * nDOF_test_element + i;
1836 for (
int j = 0; j < nDOF_trial_element; j++) {
1838 globalJacobian.data()[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += fluxJacobian_u_u[j] * u_test_dS[i];
1841 globalJacobian.data()[csrRowIndeces_w_n.data()[eN_i] + csrColumnOffsets_eb_w_n.data()[ebN_i_j]] += fluxJacobian_u_n[j] * u_test_dS[i];
1862 for (
int eN = 0; eN < nElements_global; eN++) {
1863 const int mat_eN = elementMaterialTypes.data()[eN];
1864 const double phi_eN = thetaR.data()[mat_eN] + thetaSR.data()[mat_eN];
1865 const double alpha_eN = alpha.data()[mat_eN];
1866 const double krn_end_eN = krn_end.data()[mat_eN];
1867 const double n_vg_eN =
n.data()[mat_eN];
1868 const double *KWs_eN = &KWs.data()[mat_eN *
nnz];
1871 double elementJacobian_n_n[nDOF_test_element][nDOF_trial_element];
1872 double elementJacobian_n_w[nDOF_test_element][nDOF_trial_element];
1873 for (
int i = 0; i < nDOF_test_element; i++)
1874 for (
int j = 0; j < nDOF_trial_element; j++) {
1875 elementJacobian_n_n[i][j] = 0.0;
1876 elementJacobian_n_w[i][j] = 0.0;
1878 for (
int k = 0; k < nQuadraturePoints_element; k++) {
1879 const int eN_k = eN * nQuadraturePoints_element + k;
1880 const int eN_nDOF_trial_element = eN * nDOF_trial_element;
1881 double jac[nSpace * nSpace], jacDet, jacInv[nSpace * nSpace], x_q, y_q, z_q;
1882 ck.calculateMapping_element(eN, k, mesh_dof.data(), mesh_l2g.data(),
1883 mesh_trial_ref.data(), mesh_grad_trial_ref.data(),
1884 jac, jacDet, jacInv, x_q, y_q, z_q);
1885 const double dV = std::fabs(jacDet) * dV_ref.data()[k];
1887 double u_grad_trial_qp[nDOF_trial_element * nSpace];
1888 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k * nDOF_trial_element * nSpace],
1889 jacInv, u_grad_trial_qp);
1892 ck.valFromDOF(u_dof_n.data(),
1893 &u_l2g.data()[eN_nDOF_trial_element],
1894 &u_trial_ref.data()[k * nDOF_trial_element], u_n);
1896 double u_w_qp = 0.0;
1897 ck.valFromDOF(u_dof.data(),
1898 &u_l2g.data()[eN_nDOF_trial_element],
1899 &u_trial_ref.data()[k * nDOF_trial_element], u_w_qp);
1900 double grad_u_w[nSpace];
1901 ck.gradFromDOF(u_dof.data(),
1902 &u_l2g.data()[eN_nDOF_trial_element],
1903 u_grad_trial_qp, grad_u_w);
1905 double grad_u_n[nSpace];
1906 ck.gradFromDOF(u_dof_n.data(),
1907 &u_l2g.data()[eN_nDOF_trial_element],
1908 u_grad_trial_qp, grad_u_n);
1916 const double z_cl_j = fmin(fmax(u_n, 1.0e-8), 1.0 - 1.0e-8);
1917 const double p_cl_j = fmax(u_w_qp, 1.0e2);
1920 const double S_g_j = fs_j.
S_g, Sa_j = 1.0 - S_g_j;
1922 const double N_j = fs_j.
rho_g*S_g_j + fs_j.
rho_a*Sa_j;
1927 const double dm_n_du_n = phi_eN * (dN_dz_j * z_cl_j + N_j);
1928 const double dm_n_du_w = phi_eN * (dN_dp_j * z_cl_j);
1930 const double S_wr_loc = thetaR.data()[mat_eN] / phi_eN;
1931 const double one_m_Sr_loc = 1.0 - S_wr_loc;
1932 const double Se_trap_L1965 = 1.0 - S_gr.data()[mat_eN] / one_m_Sr_loc;
1933 const double Se_raw_j = (Sa_j - S_wr_loc) / one_m_Sr_loc;
1934 double Se_a, dSe_dp, dSe_dz;
1935 if (Se_raw_j <= 0.0) { Se_a = 0.0; dSe_dp = 0.0; dSe_dz = 0.0; }
1936 else if (Se_raw_j >= 1.0) { Se_a = 1.0; dSe_dp = 0.0; dSe_dz = 0.0; }
1937 else { Se_a = Se_raw_j; dSe_dp = -fs_j.
dS_g_dp/one_m_Sr_loc; dSe_dz = -fs_j.
dS_g_dz/one_m_Sr_loc; }
1938 double KWr_a=0.0, DKWr_a=0.0, thW_a=0.0, DthW_a=0.0, KNr_a=0.0, DKNr_a=0.0;
1939 double pc_a=0.0, dpc_dSe_a=0.0, d2pc_a=0.0;
1942 thetaR.data()[mat_eN], thetaSR.data()[mat_eN], thW_a, DthW_a, KWr_a, DKWr_a);
1947 thetaR.data()[mat_eN], thetaSR.data()[mat_eN], thW_a, DthW_a, KWr_a, DKWr_a);
1951 KNr_a *= krn_end_eN; DKNr_a *= krn_end_eN;
1952 const double pcp_a = dpc_dSe_a / one_m_Sr_loc;
1953 const double dpcp_dp = (d2pc_a / one_m_Sr_loc) * dSe_dp;
1954 const double dpcp_dz = (d2pc_a / one_m_Sr_loc) * dSe_dz;
1959 const double rho_g_mass = fs_j.
rho_g*Mbar_g, rho_a_mass = fs_j.
rho_a*Mbar_a;
1961 const double drgm_dz = fs_j.
rho_g*fs_j.
dY_dz*dMm;
1965 const double Ag = fs_j.
rho_g*fs_j.
Y, Aa = fs_j.
rho_a*fs_j.
X;
1972 double ug[nSpace], ua[nSpace];
1973 double dug_dp[nSpace], dug_dz[nSpace], dua_dp[nSpace], dua_dz[nSpace];
1974 for (
int I = 0; I < nSpace; I++) {
1975 double ugI=0.0, uaI=0.0, dugp=0.0, dugz=0.0, duap=0.0, duaz=0.0;
1976 for (
int ii = a_rowptr.data()[I]; ii < a_rowptr.data()[I + 1]; ii++) {
1977 const int J = a_colind.data()[ii];
1978 const double Kii = KWs_eN[ii];
1979 const double Mob_g = KNr_a*Kii/mu_n, Mob_a = KWr_a*Kii;
1980 const double dMobg_dp = (DKNr_a*Kii/mu_n)*dSe_dp, dMobg_dz = (DKNr_a*Kii/mu_n)*dSe_dz;
1981 const double dMoba_dp = (DKWr_a*Kii)*dSe_dp, dMoba_dz = (DKWr_a*Kii)*dSe_dz;
1982 const double gJ = gravity.data()[J];
1983 const double gradSa = -(fs_j.
dS_g_dp*grad_u_w[J] + fs_j.
dS_g_dz*grad_u_n[J]);
1984 const double gp_a = grad_u_w[J] - rho_a_mass*gJ;
1985 const double gp_g = grad_u_w[J] + pcp_a*gradSa - rho_g_mass*gJ;
1986 ugI -= Mob_g*gp_g; uaI -= Mob_a*gp_a;
1989 const double dgpg_dp = dpcp_dp*gradSa + pcp_a*dgradSa_dp - drgm_dp*gJ;
1990 const double dgpg_dz = dpcp_dz*gradSa + pcp_a*dgradSa_dz - drgm_dz*gJ;
1991 dugp -= dMobg_dp*gp_g + Mob_g*dgpg_dp;
1992 dugz -= dMobg_dz*gp_g + Mob_g*dgpg_dz;
1993 duap -= dMoba_dp*gp_a + Mob_a*(-dram_dp*gJ);
1994 duaz -= dMoba_dz*gp_a + Mob_a*(-dram_dz*gJ);
1996 ug[I]=ugI; ua[I]=uaI;
1997 dug_dp[I]=dugp; dug_dz[I]=dugz; dua_dp[I]=duap; dua_dz[I]=duaz;
2005 for (
int i = 0; i < nDOF_test_element; i++) {
2006 const double test_i = u_test_ref.data()[k * nDOF_test_element + i];
2008 double Sval_p = 0.0, Sval_z = 0.0;
2009 for (
int I = 0; I < nSpace; I++) {
2010 const double gNiI = u_grad_trial_qp[i * nSpace + I];
2011 Sval_p += (dAg_dp*ug[I] + Ag*dug_dp[I] + dAa_dp*ua[I] + Aa*dua_dp[I]) * gNiI;
2012 Sval_z += (dAg_dz*ug[I] + Ag*dug_dz[I] + dAa_dz*ua[I] + Aa*dua_dz[I]) * gNiI;
2014 for (
int j = 0; j < nDOF_trial_element; j++) {
2015 const double trial_j = u_trial_ref.data()[k * nDOF_trial_element + j];
2017 double Sgrad_p = 0.0, Sgrad_z = 0.0;
2018 for (
int I = 0; I < nSpace; I++) {
2019 const double gNiI = u_grad_trial_qp[i * nSpace + I];
2020 for (
int ii = a_rowptr.data()[I]; ii < a_rowptr.data()[I + 1]; ii++) {
2021 const int J = a_colind.data()[ii];
2022 const double Kii = KWs_eN[ii];
2023 const double Mob_g = KNr_a*Kii/mu_n, Mob_a = KWr_a*Kii;
2024 const double gNjJ = u_grad_trial_qp[j * nSpace + J];
2026 const double dFdgp = Ag*(-Mob_g*(1.0 - pcp_a*fs_j.
dS_g_dp)) + Aa*(-Mob_a);
2028 const double dFdgz = Ag*(-Mob_g*(-pcp_a*fs_j.
dS_g_dz));
2029 Sgrad_p += dFdgp * gNjJ * gNiI;
2030 Sgrad_z += dFdgz * gNjJ * gNiI;
2034 elementJacobian_n_n[i][j] += (dm_n_du_n * test_i * trial_j * dV) / dt_n;
2035 elementJacobian_n_n[i][j] -= (Sval_z * trial_j + Sgrad_z) * dV;
2037 elementJacobian_n_w[i][j] += (dm_n_du_w * test_i * trial_j * dV) / dt_n;
2038 elementJacobian_n_w[i][j] -= (Sval_p * trial_j + Sgrad_p) * dV;
2042 for (
int i = 0; i < nDOF_test_element; i++) {
2043 const int eN_i = eN * nDOF_test_element + i;
2044 for (
int j = 0; j < nDOF_trial_element; j++) {
2045 const int eN_i_j = eN_i * nDOF_trial_element + j;
2046 globalJacobian.data()[csrRowIndeces_n_n.data()[eN_i] + csrColumnOffsets_n_n.data()[eN_i_j]]
2047 += elementJacobian_n_n[i][j];
2048 globalJacobian.data()[csrRowIndeces_n_w.data()[eN_i] + csrColumnOffsets_n_w.data()[eN_i_j]]
2049 += elementJacobian_n_w[i][j];
2065 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++) {
2066 const int ebN = exteriorElementBoundariesArray.data()[ebNE];
2067 const int eN = elementBoundaryElementsArray.data()[ebN * 2 + 0];
2068 const int ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN * 2 + 0];
2069 const int eN_nDOF_trial_element = eN * nDOF_trial_element;
2070 const int mat_eN = elementMaterialTypes.data()[eN];
2071 const double phi_eN = thetaR.data()[mat_eN] + thetaSR.data()[mat_eN];
2072 const double alpha_eN = alpha.data()[mat_eN];
2073 const double krn_end_eN = krn_end.data()[mat_eN];
2074 const double n_vg_eN =
n.data()[mat_eN];
2075 const double *KWs_eN = &KWs.data()[mat_eN *
nnz];
2076 const double S_wr_loc = thetaR.data()[mat_eN] / phi_eN;
2077 const double one_m_Sr_loc = 1.0 - S_wr_loc;
2078 const double Se_trap_L2110 = 1.0 - S_gr.data()[mat_eN] / one_m_Sr_loc;
2080 double elementJacobian_n_n_eb[nDOF_test_element][nDOF_trial_element];
2081 double elementJacobian_n_w_eb[nDOF_test_element][nDOF_trial_element];
2082 for (
int i = 0; i < nDOF_test_element; i++)
2083 for (
int j = 0; j < nDOF_trial_element; j++) {
2084 elementJacobian_n_n_eb[i][j] = 0.0;
2085 elementJacobian_n_w_eb[i][j] = 0.0;
2088 for (
int kb = 0; kb < nQuadraturePoints_elementBoundary; kb++) {
2089 const int ebNE_kb = ebNE * nQuadraturePoints_elementBoundary + kb;
2090 const int ebN_local_kb = ebN_local * nQuadraturePoints_elementBoundary + kb;
2091 const int ebN_local_kb_nSpace = ebN_local_kb * nSpace;
2093 double jac_ext[nSpace * nSpace], jacDet_ext, jacInv_ext[nSpace * nSpace];
2094 double boundaryJac_b[nSpace * (nSpace - 1)];
2095 double metricTensor_b[(nSpace - 1) * (nSpace - 1)];
2096 double metricTensorDetSqrt_b, dS_eb, normal_b[3];
2097 double xt_b, yt_b, zt_b, integralScaling_b;
2098 double x_eb, y_eb, z_eb;
2099 ck.calculateMapping_elementBoundary(eN, ebN_local, kb, ebN_local_kb,
2100 mesh_dof.data(), mesh_l2g.data(), mesh_trial_trace_ref.data(),
2101 mesh_grad_trial_trace_ref.data(), boundaryJac_ref.data(),
2102 jac_ext, jacDet_ext, jacInv_ext, boundaryJac_b, metricTensor_b,
2103 metricTensorDetSqrt_b, normal_ref.data(), normal_b,
2105 ck.calculateMappingVelocity_elementBoundary(eN, ebN_local, kb, ebN_local_kb,
2106 mesh_velocity_dof.data(), mesh_l2g.data(), mesh_trial_trace_ref.data(),
2107 xt_b, yt_b, zt_b, normal_b, boundaryJac_b, metricTensor_b,
2109 dS_eb = ((1.0 - MOVING_DOMAIN) * metricTensorDetSqrt_b
2110 + MOVING_DOMAIN * integralScaling_b) * dS_ref.data()[kb];
2112 double u_grad_trial_trace_b[nDOF_trial_element * nSpace];
2113 ck.gradTrialFromRef(
2114 &u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace * nDOF_trial_element],
2115 jacInv_ext, u_grad_trial_trace_b);
2116 double u_w_ext_b = 0.0, u_n_ext_b = 0.0;
2117 double grad_u_w_ext_b[nSpace], grad_u_n_ext_b[nSpace];
2118 ck.valFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element],
2119 &u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element], u_w_ext_b);
2120 ck.valFromDOF(u_dof_n.data(), &u_l2g.data()[eN_nDOF_trial_element],
2121 &u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element], u_n_ext_b);
2122 ck.gradFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element],
2123 u_grad_trial_trace_b, grad_u_w_ext_b);
2124 ck.gradFromDOF(u_dof_n.data(), &u_l2g.data()[eN_nDOF_trial_element],
2125 u_grad_trial_trace_b, grad_u_n_ext_b);
2127 const int isDir_n = isDOFBoundary_n.data()[ebNE_kb];
2128 const double penalty = ebqe_penalty_ext.data()[ebNE_kb];
2138 const double z_clb = fmin(fmax(u_n_ext_b, 1.0e-8), 1.0 - 1.0e-8);
2139 const double p_clb = fmax(u_w_ext_b, 1.0e2);
2142 const double Sa_b = 1.0 - fsb.
S_g;
2143 const double Se_raw_b = (Sa_b - S_wr_loc) / one_m_Sr_loc;
2144 double Se_b, dSe_dp_b, dSe_dz_b;
2145 if (Se_raw_b <= 0.0) { Se_b = 0.0; dSe_dp_b = 0.0; dSe_dz_b = 0.0; }
2146 else if (Se_raw_b >= 1.0) { Se_b = 1.0; dSe_dp_b = 0.0; dSe_dz_b = 0.0; }
2147 else { Se_b = Se_raw_b; dSe_dp_b = -fsb.
dS_g_dp/one_m_Sr_loc; dSe_dz_b = -fsb.
dS_g_dz/one_m_Sr_loc; }
2148 double KWr_b=0,DKWr_b=0,thW_b=0,DthW_b=0,KNr_b=0,DKNr_b=0,pc_b=0,dpc_dSe_b=0,d2pc_b=0;
2158 KNr_b *= krn_end_eN; DKNr_b *= krn_end_eN;
2159 const double pcp_b = dpc_dSe_b / one_m_Sr_loc;
2160 const double dpcp_dp_b = (d2pc_b / one_m_Sr_loc) * dSe_dp_b;
2161 const double dpcp_dz_b = (d2pc_b / one_m_Sr_loc) * dSe_dz_b;
2165 const double rho_g_mass_b = fsb.
rho_g*Mbar_g_b, rho_a_mass_b = fsb.
rho_a*Mbar_a_b;
2167 const double drgm_dz_b = fsb.
rho_g*fsb.
dY_dz*dMm_b;
2170 const double Ag = fsb.
rho_g*fsb.
Y, Aa = fsb.
rho_a*fsb.
X;
2175 double ug_b[nSpace], ua_b[nSpace];
2176 double dug_dp_b[nSpace], dug_dz_b[nSpace], dua_dp_b[nSpace], dua_dz_b[nSpace];
2177 for (
int I = 0; I < nSpace; I++) {
2178 double ugI=0.0, uaI=0.0, dugp=0.0, dugz=0.0, duap=0.0, duaz=0.0;
2179 for (
int ii = a_rowptr.data()[I]; ii < a_rowptr.data()[I + 1]; ii++) {
2180 const int J = a_colind.data()[ii];
2181 const double Kii = KWs_eN[ii];
2182 const double Mob_g = KNr_b*Kii/mu_n, Mob_a = KWr_b*Kii;
2183 const double dMobg_dp = (DKNr_b*Kii/mu_n)*dSe_dp_b, dMobg_dz = (DKNr_b*Kii/mu_n)*dSe_dz_b;
2184 const double dMoba_dp = (DKWr_b*Kii)*dSe_dp_b, dMoba_dz = (DKWr_b*Kii)*dSe_dz_b;
2185 const double gJ = gravity.data()[J];
2186 const double gradSa = -(fsb.
dS_g_dp*grad_u_w_ext_b[J] + fsb.
dS_g_dz*grad_u_n_ext_b[J]);
2187 const double gp_a = grad_u_w_ext_b[J] - rho_a_mass_b*gJ;
2188 const double gp_g = grad_u_w_ext_b[J] + pcp_b*gradSa - rho_g_mass_b*gJ;
2189 ugI -= Mob_g*gp_g; uaI -= Mob_a*gp_a;
2190 const double dgradSa_dp = -(fsb.
d2S_g_dp2 *grad_u_w_ext_b[J] + fsb.
d2S_g_dpdz*grad_u_n_ext_b[J]);
2191 const double dgradSa_dz = -(fsb.
d2S_g_dpdz*grad_u_w_ext_b[J] + fsb.
d2S_g_dz2 *grad_u_n_ext_b[J]);
2192 const double dgpg_dp = dpcp_dp_b*gradSa + pcp_b*dgradSa_dp - drgm_dp_b*gJ;
2193 const double dgpg_dz = dpcp_dz_b*gradSa + pcp_b*dgradSa_dz - drgm_dz_b*gJ;
2194 dugp -= dMobg_dp*gp_g + Mob_g*dgpg_dp;
2195 dugz -= dMobg_dz*gp_g + Mob_g*dgpg_dz;
2196 duap -= dMoba_dp*gp_a + Mob_a*(-dram_dp_b*gJ);
2197 duaz -= dMoba_dz*gp_a + Mob_a*(-dram_dz_b*gJ);
2199 ug_b[I]=ugI; ua_b[I]=uaI;
2200 dug_dp_b[I]=dugp; dug_dz_b[I]=dugz; dua_dp_b[I]=duap; dua_dz_b[I]=duaz;
2203 double Sval_p_b = 0.0, Sval_z_b = 0.0;
2204 for (
int I = 0; I < nSpace; I++) {
2205 Sval_p_b += (dAg_dp*ug_b[I] + Ag*dug_dp_b[I] + dAa_dp*ua_b[I] + Aa*dua_dp_b[I]) * normal_b[I];
2206 Sval_z_b += (dAg_dz*ug_b[I] + Ag*dug_dz_b[I] + dAa_dz*ua_b[I] + Aa*dua_dz_b[I]) * normal_b[I];
2210 double Kw_rep = 0.0;
2211 for (
int ii = 0; ii <
nnz; ii++) Kw_rep = fmax(Kw_rep, fabs(KWs_eN[ii]));
2212 const double a_n_scale = (Ag*KNr_b/mu_n + Aa*KWr_b) * Kw_rep;
2213 const double pen_eff = penalty * a_n_scale;
2215 for (
int i = 0; i < nDOF_test_element; i++) {
2216 const double test_i_dS = u_test_trace_ref.data()[
2217 ebN_local_kb * nDOF_test_element + i] * dS_eb;
2218 for (
int j = 0; j < nDOF_trial_element; j++) {
2219 const double trial_j_b = u_trial_trace_ref.data()[
2220 ebN_local_kb * nDOF_test_element + j];
2221 double Sgrad_p_b = 0.0, Sgrad_z_b = 0.0;
2222 for (
int I = 0; I < nSpace; I++) {
2223 for (
int ii = a_rowptr.data()[I]; ii < a_rowptr.data()[I + 1]; ii++) {
2224 const int J = a_colind.data()[ii];
2225 const double Kii = KWs_eN[ii];
2226 const double Mob_g = KNr_b*Kii/mu_n, Mob_a = KWr_b*Kii;
2227 const double gNjJ = u_grad_trial_trace_b[j * nSpace + J];
2228 const double dFdgp = Ag*(-Mob_g*(1.0 - pcp_b*fsb.
dS_g_dp)) + Aa*(-Mob_a);
2229 const double dFdgz = Ag*(-Mob_g*(-pcp_b*fsb.
dS_g_dz));
2230 Sgrad_p_b += dFdgp * gNjJ * normal_b[I];
2231 Sgrad_z_b += dFdgz * gNjJ * normal_b[I];
2234 double jac_nn = Sval_z_b * trial_j_b + Sgrad_z_b;
2235 double jac_nw = Sval_p_b * trial_j_b + Sgrad_p_b;
2237 jac_nn += pen_eff * trial_j_b;
2242 elementJacobian_n_n_eb[i][j] += jac_nn * test_i_dS;
2243 elementJacobian_n_w_eb[i][j] += jac_nw * test_i_dS;
2248 for (
int i = 0; i < nDOF_test_element; i++) {
2249 const int eN_i = eN * nDOF_test_element + i;
2250 for (
int j = 0; j < nDOF_trial_element; j++) {
2252 + i * nDOF_trial_element + j;
2253 globalJacobian.data()[csrRowIndeces_n_n.data()[eN_i]
2254 + csrColumnOffsets_eb_n_n.data()[ebN_i_j]] += elementJacobian_n_n_eb[i][j];
2255 globalJacobian.data()[csrRowIndeces_n_w.data()[eN_i]
2256 + csrColumnOffsets_eb_n_w.data()[ebN_i_j]] += elementJacobian_n_w_eb[i][j];
2632 xt::pyarray<double> &globalJacobian = args.
array<
double>(
"globalJacobian");
2633 double Theta = args.
scalar<
double>(
"Theta");
2634 double Theta_h = args.
scalar<
double>(
"Theta_h");
2635 xt::pyarray<double> &bc_mask = args.
array<
double>(
"bc_mask");
2636 double dt = args.
scalar<
double>(
"dt");
2637 xt::pyarray<double> &mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
2638 xt::pyarray<double> &mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
2639 xt::pyarray<double> &mesh_dof = args.
array<
double>(
"mesh_dof");
2640 xt::pyarray<double> &mesh_velocity_dof = args.
array<
double>(
"mesh_velocity_dof");
2641 double MOVING_DOMAIN = args.
scalar<
double>(
"MOVING_DOMAIN");
2642 xt::pyarray<int> &mesh_l2g = args.
array<
int>(
"mesh_l2g");
2643 xt::pyarray<double> &dV_ref = args.
array<
double>(
"dV_ref");
2644 xt::pyarray<double> &u_trial_ref = args.
array<
double>(
"u_trial_ref");
2645 xt::pyarray<double> &u_grad_trial_ref = args.
array<
double>(
"u_grad_trial_ref");
2646 xt::pyarray<double> &u_test_ref = args.
array<
double>(
"u_test_ref");
2647 xt::pyarray<double> &u_grad_test_ref = args.
array<
double>(
"u_grad_test_ref");
2648 xt::pyarray<double> &mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
2649 xt::pyarray<double> &mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
2650 xt::pyarray<double> &dS_ref = args.
array<
double>(
"dS_ref");
2651 xt::pyarray<double> &u_trial_trace_ref = args.
array<
double>(
"u_trial_trace_ref");
2653 xt::pyarray<double> &u_grad_trial_trace_ref = args.
array<
double>(
"u_grad_trial_trace_ref");
2654 xt::pyarray<double> &u_test_trace_ref = args.
array<
double>(
"u_test_trace_ref");
2655 xt::pyarray<double> &u_grad_test_trace_ref = args.
array<
double>(
"u_grad_test_trace_ref");
2656 xt::pyarray<double> &normal_ref = args.
array<
double>(
"normal_ref");
2657 xt::pyarray<double> &boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
2658 int nElements_global = args.
scalar<
int>(
"nElements_global");
2659 xt::pyarray<double> &ebqe_penalty_ext = args.
array<
double>(
"ebqe_penalty_ext");
2660 xt::pyarray<int> &elementMaterialTypes = args.
array<
int>(
"elementMaterialTypes");
2661 xt::pyarray<int> &isSeepageFace = args.
array<
int>(
"isSeepageFace");
2662 xt::pyarray<int> &a_rowptr = args.
array<
int>(
"a_rowptr");
2663 xt::pyarray<int> &a_colind = args.
array<
int>(
"a_colind");
2664 double rho = args.
scalar<
double>(
"rho");
2665 double beta = args.
scalar<
double>(
"beta");
2667 xt::pyarray<double> &q_rho = args.
array<
double>(
"q_rho");
2668 xt::pyarray<double> &ebqe_rho = args.
array<
double>(
"ebqe_rho");
2672 xt::pyarray<double> &gravity = args.
array<
double>(
"gravity");
2673 xt::pyarray<double> &alpha = args.
array<
double>(
"alpha");
2674 xt::pyarray<double> &
n = args.
array<
double>(
"n");
2675 xt::pyarray<double> &thetaR = args.
array<
double>(
"thetaR");
2676 xt::pyarray<double> &thetaSR = args.
array<
double>(
"thetaSR");
2677 xt::pyarray<double> &KWs = args.
array<
double>(
"KWs");
2678 xt::pyarray<double> &krn_end = args.
array<
double>(
"krn_end");
2679 xt::pyarray<double> &S_gr = args.
array<
double>(
"S_gr");
2680 double mu_n = args.
scalar<
double>(
"mu_n");
2681 double useMetrics = args.
scalar<
double>(
"useMetrics");
2682 double alphaBDF = args.
scalar<
double>(
"alphaBDF");
2683 int lag_shockCapturing = args.
scalar<
int>(
"lag_shockCapturing");
2684 double shockCapturingDiffusion = args.
scalar<
double>(
"shockCapturingDiffusion");
2685 double sc_uref = args.
scalar<
double>(
"sc_uref");
2686 double sc_alpha = args.
scalar<
double>(
"sc_alpha");
2687 xt::pyarray<int> &u_l2g = args.
array<
int>(
"u_l2g");
2697 xt::pyarray<int> &u_l2g_n = args.
array<
int>(
"u_l2g_n");
2698 const int split_z = args.
scalar<
int>(
"split_z");
2699 const double D_m = args.
scalar<
double>(
"D_m");
2700 xt::pyarray<int> &interface_pairs = args.
array<
int>(
"interface_pairs");
2701 const int n_interface_pairs = args.
scalar<
int>(
"n_interface_pairs");
2706 const double split_anchor_alpha = args.
scalar<
double>(
"split_anchor_alpha");
2711 const double split_anchor_Sg_tol = args.
scalar<
double>(
"split_anchor_Sg_tol");
2712 const double split_anchor_X_tol = args.
scalar<
double>(
"split_anchor_X_tol");
2714 const double split_anchor_zfloor = args.
scalar<
double>(
"split_anchor_zfloor");
2715 (void)split_anchor_zfloor;
2721 const int split_anchor_layer1 = args.
scalar<
int>(
"split_anchor_layer1");
2722 xt::pyarray<int> &r_l2g = args.
array<
int>(
"r_l2g");
2723 xt::pyarray<double> &elementDiameter = args.
array<
double>(
"elementDiameter");
2724 int degree_polynomial = args.
scalar<
int>(
"degree_polynomial");
2725 xt::pyarray<double> &u_dof = args.
array<
double>(
"u_dof");
2726 xt::pyarray<double> &u_dof_old = args.
array<
double>(
"u_dof_old");
2727 xt::pyarray<double> &velocity = args.
array<
double>(
"velocity");
2728 xt::pyarray<double> &q_m = args.
array<
double>(
"q_m");
2729 xt::pyarray<double> &q_theta = args.
array<
double>(
"q_theta");
2730 xt::pyarray<double> &q_u = args.
array<
double>(
"q_u");
2731 xt::pyarray<double> &q_dV = args.
array<
double>(
"q_dV");
2732 xt::pyarray<double> &q_m_betaBDF = args.
array<
double>(
"q_m_betaBDF");
2733 xt::pyarray<double> &cfl = args.
array<
double>(
"cfl");
2734 xt::pyarray<double> &q_numDiff_u = args.
array<
double>(
"q_numDiff_u");
2735 xt::pyarray<double> &q_numDiff_u_last = args.
array<
double>(
"q_numDiff_u_last");
2736 int offset_u = args.
scalar<
int>(
"offset_u");
2737 int stride_u = args.
scalar<
int>(
"stride_u");
2740 xt::pyarray<double> &u_dof_n = args.
array<
double>(
"u_dof_n");
2741 xt::pyarray<double> &u_dof_n_old = args.
array<
double>(
"u_dof_n_old");
2752 const double rho_n = args.
scalar<
double>(
"rho_n");
2753 const double p_ref_n = args.
scalar<
double>(
"p_ref_n");
2754 const bool rho_n_compressible = (p_ref_n > 0.0);
2755 const double inv_p_ref_n = rho_n_compressible ? (1.0 / p_ref_n) : 0.0;
2756 const int offset_n = args.
scalar<
int>(
"offset_n");
2757 const int stride_n = args.
scalar<
int>(
"stride_n");
2764 const int inj_point_mode = args.
scalar<
int>(
"inj_point_mode");
2765 const int inj_n_ports = args.
scalar<
int>(
"inj_n_ports");
2766 xt::pyarray<int> &inj_element = args.
array<
int>(
"inj_element");
2767 xt::pyarray<double> &inj_weight = args.
array<
double>(
"inj_weight");
2768 xt::pyarray<double> &inj_rate = args.
array<
double>(
"inj_rate");
2774 xt::pyarray<double> &c_dof = args.
array<
double>(
"c_dof");
2775 const double k_d = args.
scalar<
double>(
"k_d");
2776 const double c_sat = args.
scalar<
double>(
"c_sat");
2780 xt::pyarray<double> &injection_dof = args.
array<
double>(
"injection_dof");
2781 xt::pyarray<double> &globalResidual = args.
array<
double>(
"globalResidual");
2782 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
2783 xt::pyarray<int> &exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
2784 xt::pyarray<int> &elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
2785 xt::pyarray<int> &elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
2786 xt::pyarray<double> &ebqe_velocity_ext = args.
array<
double>(
"ebqe_velocity_ext");
2787 xt::pyarray<int> &isDOFBoundary_u = args.
array<
int>(
"isDOFBoundary_u");
2788 xt::pyarray<double> &ebqe_bc_u_ext = args.
array<
double>(
"ebqe_bc_u_ext");
2790 xt::pyarray<int> &isDOFBoundary_n = args.
array<
int>(
"isDOFBoundary_n");
2791 xt::pyarray<double> &ebqe_bc_u_n_ext = args.
array<
double>(
"ebqe_bc_u_n_ext");
2792 xt::pyarray<int> &isFluxBoundary_u = args.
array<
int>(
"isFluxBoundary_u");
2793 xt::pyarray<double> &ebqe_bc_flux_ext = args.
array<
double>(
"ebqe_bc_flux_ext");
2794 xt::pyarray<double> &ebqe_phi = args.
array<
double>(
"ebqe_phi");
2795 double epsFact = args.
scalar<
double>(
"epsFact");
2796 xt::pyarray<double> &ebqe_u = args.
array<
double>(
"ebqe_u");
2797 xt::pyarray<double> &ebqe_theta = args.
array<
double>(
"ebqe_theta");
2798 xt::pyarray<double> &ebqe_flux = args.
array<
double>(
"ebqe_flux");
2800 double cE = args.
scalar<
double>(
"cE");
2801 double cK = args.
scalar<
double>(
"cK");
2803 double uL = args.
scalar<
double>(
"uL");
2804 double uR = args.
scalar<
double>(
"uR");
2806 int numDOFs = args.
scalar<
int>(
"numDOFs");
2810 int numDOFs_u = args.
scalar<
int>(
"numDOFs_u");
2811 int NNZ = args.
scalar<
int>(
"NNZ");
2812 xt::pyarray<int> &csrRowIndeces_DofLoops = args.
array<
int>(
"csrRowIndeces_DofLoops");
2813 xt::pyarray<int> &csrColumnOffsets_DofLoops = args.
array<
int>(
"csrColumnOffsets_DofLoops");
2814 xt::pyarray<int> &csrRowIndeces_Full = args.
array<
int>(
"csrRowIndeces_Full");
2815 xt::pyarray<int> &csrColumnOffsets_Full = args.
array<
int>(
"csrColumnOffsets_Full");
2816 xt::pyarray<int> &csrRowIndeces_CellLoops = args.
array<
int>(
"csrRowIndeces_CellLoops");
2817 xt::pyarray<int> &csrColumnOffsets_CellLoops = args.
array<
int>(
"csrColumnOffsets_CellLoops");
2818 xt::pyarray<int> &csrColumnOffsets_eb_CellLoops = args.
array<
int>(
"csrColumnOffsets_eb_CellLoops");
2820 xt::pyarray<double> &Cx = args.
array<
double>(
"Cx");
2821 xt::pyarray<double> &Cy = args.
array<
double>(
"Cy");
2822 xt::pyarray<double> &Cz = args.
array<
double>(
"Cz");
2823 xt::pyarray<double> &CTx = args.
array<
double>(
"CTx");
2824 xt::pyarray<double> &CTy = args.
array<
double>(
"CTy");
2825 xt::pyarray<double> &CTz = args.
array<
double>(
"CTz");
2826 xt::pyarray<double> &ML = args.
array<
double>(
"ML");
2827 xt::pyarray<double> &MC = args.
array<
double>(
"MC");
2829 xt::pyarray<double> &delta_x_ij = args.
array<
double>(
"delta_x_ij");
2831 int LUMPED_MASS_MATRIX = args.
scalar<
int>(
"LUMPED_MASS_MATRIX");
2834 int ENTROPY_TYPE = args.
scalar<
int>(
"ENTROPY_TYPE");
2839 xt::pyarray<double> &dLow = args.
array<
double>(
"dLow");
2840 xt::pyarray<double> &fluxMatrix = args.
array<
double>(
"fluxMatrix");
2841 xt::pyarray<double> &mDotLow = args.
array<
double>(
"mDotLow");
2842 xt::pyarray<double> &mLow = args.
array<
double>(
"mLow");
2843 xt::pyarray<double> &dt_times_fH_minus_fL = args.
array<
double>(
"dt_times_fH_minus_fL");
2844 xt::pyarray<double> &min_m_bc = args.
array<
double>(
"min_m_bc");
2845 xt::pyarray<double> &max_m_bc = args.
array<
double>(
"max_m_bc");
2847 xt::pyarray<double> &quantDOFs = args.
array<
double>(
"quantDOFs");
2848 xt::pyarray<double> &mn = args.
array<
double>(
"mn");
2849 xt::pyarray<double> &fluxCorrection = args.
array<
double>(
"fluxCorrection");
2850 xt::pyarray<double> &limited_solution = args.
array<
double>(
"limited_solution");
2851 xt::pyarray<int> &freeDOFMaterialTypes = args.
array<
int>(
"freeDOFMaterialTypes");
2852 xt::pyarray<int> &freeDOFToNode_u = args.
array<
int>(
"freeDOFToNode_u");
2859 xt::pyarray<int> &node2zdof = args.
array<
int>(
"node2zdof");
2865 xt::pyarray<int> &nodeMaterialTypes_n = args.
array<
int>(
"nodeMaterialTypes_n");
2868 xt::pyarray<double> &node_pd_min = args.
array<
double>(
"node_pd_min");
2871 xt::pyarray<double> &node_Sn_max = args.
array<
double>(
"node_Sn_max");
2875 xt::pyarray<double> &gas_diag = args.
array<
double>(
"gas_diag");
2877 xt::pyarray<double> &velocity_couple = args.
array<
double>(
"velocity_couple");
2878 xt::pyarray<double> &ebqe_velocity_ext_couple = args.
array<
double>(
"ebqe_velocity_ext_couple");
2882 xt::pyarray<double> &anb_seepage_flux_n = args.
array<
double>(
"anb_seepage_flux_n");
2883 xt::pyarray<double> &q_velocity = args.
array<
double>(
"q_velocity");
2884 double &anb_seepage_flux(args.
scalar<
double>(
"anb_seepage_flux"));
2885 anb_seepage_flux = 0.0;
2886 xt::pyarray<int> &csrRowIndeces_u_u = args.
array<
int>(
"csrRowIndeces_u_u");
2887 xt::pyarray<int> &csrColumnOffsets_u_u = args.
array<
int>(
"csrColumnOffsets_u_u");
2888 xt::pyarray<int> &csrColumnOffsets_eb_u_u = args.
array<
int>(
"csrColumnOffsets_eb_u_u");
2892 xt::pyarray<int> &csrRowIndeces_n_n = args.
array<
int>(
"csrRowIndeces_n_n");
2896 xt::pyarray<int> &csrRowIndeces_n_w = args.
array<
int>(
"csrRowIndeces_n_w");
2897 xt::pyarray<int> &csrColumnOffsets_n_n = args.
array<
int>(
"csrColumnOffsets_n_n");
2898 xt::pyarray<int> &csrColumnOffsets_n_w = args.
array<
int>(
"csrColumnOffsets_n_w");
2904 xt::pyarray<int> &csrRowIndeces_w_w = args.
array<
int>(
"csrRowIndeces_w_w");
2905 xt::pyarray<int> &csrColumnOffsets_w_w = args.
array<
int>(
"csrColumnOffsets_w_w");
2906 xt::pyarray<int> &csrRowIndeces_w_n = args.
array<
int>(
"csrRowIndeces_w_n");
2907 xt::pyarray<int> &csrColumnOffsets_w_n = args.
array<
int>(
"csrColumnOffsets_w_n");
2910 xt::pyarray<int> &csrColumnOffsets_eb_n_n = args.
array<
int>(
"csrColumnOffsets_eb_n_n");
2911 xt::pyarray<int> &csrColumnOffsets_eb_n_w = args.
array<
int>(
"csrColumnOffsets_eb_n_w");
2917 int numDOFs_n = args.
scalar<
int>(
"numDOFs_n");
2918 int NNZ_n = args.
scalar<
int>(
"NNZ_n");
2919 xt::pyarray<int> &csrRowIndeces_n_DofLoops = args.
array<
int>(
"csrRowIndeces_n_DofLoops");
2920 xt::pyarray<int> &csrColumnOffsets_n_DofLoops = args.
array<
int>(
"csrColumnOffsets_n_DofLoops");
2924 xt::pyarray<int> &comp1_full_offsets = args.
array<
int>(
"comp1_full_offsets");
2932 xt::pyarray<int> &comp10_full_offsets = args.
array<
int>(
"comp10_full_offsets");
2939 xt::pyarray<int> &comp1_iface_offsets = args.
array<
int>(
"comp1_iface_offsets");
2940 xt::pyarray<double> &dLow_n = args.
array<
double>(
"dLow_n");
2941 xt::pyarray<double> &dEV_n = args.
array<
double>(
"dEV_n");
2942 xt::pyarray<double> &fluxMatrix_n = args.
array<
double>(
"fluxMatrix_n");
2943 xt::pyarray<double> &mLow_n = args.
array<
double>(
"mLow_n");
2944 xt::pyarray<double> &mDotLow_n = args.
array<
double>(
"mDotLow_n");
2945 double u_n_L = args.
scalar<
double>(
"u_n_L");
2946 double u_n_R = args.
scalar<
double>(
"u_n_R");
2947 xt::pyarray<double> &mn_n = args.
array<
double>(
"mn_n");
2948 xt::pyarray<double> &quantDOFs_n = args.
array<
double>(
"quantDOFs_n");
2962 xt::pyarray<double> &gas_budget_node = args.
array<
double>(
"gas_budget_node");
2965 xt::pyarray<double> &dt_times_fH_minus_fL_n = args.
array<
double>(
"dt_times_fH_minus_fL_n");
2966 xt::pyarray<double> &fluxCorrection_n = args.
array<
double>(
"fluxCorrection_n");
2967 int FCT_n = args.
scalar<
int>(
"FCT_n");
2969 std::vector<double> Rpos(numDOFs, 0.0), Rneg(numDOFs, 0.0);
2970 std::vector<double> TransportMatrix(NNZ, 0.0),
2971 TransportMatrixConsistent(NNZ, 0.0),
2972 TransportMatrixn(NNZ, 0.0),
2973 TransportMatrixConsistentn(NNZ, 0.0);
2980 std::valarray<double> u_free_dof(numDOFs);
2981 std::valarray<double> u_free_dof_old(numDOFs);
2982 std::valarray<double> ML2(numDOFs);
2984 std::vector<double> rho_dof(numDOFs, 0.0);
2985 std::vector<double> ML_rho(numDOFs, 0.0);
2986 std::fill(velocity_couple.data(), velocity_couple.data() + velocity_couple.size(), 0.0);
2987 std::fill(ebqe_velocity_ext_couple.data(), ebqe_velocity_ext_couple.data() + ebqe_velocity_ext_couple.size(), 0.0);
2988 auto full_offset_from_compact = [&](
int i_compact,
int j_compact) ->
int
2990 const int full_i = offset_u + stride_u * i_compact;
2991 const int full_j = offset_u + stride_u * j_compact;
2992 for (
int offset = csrRowIndeces_Full.data()[full_i]; offset < csrRowIndeces_Full.data()[full_i + 1]; ++offset)
2993 if (csrColumnOffsets_Full.data()[offset] == full_j)
return offset;
2997 for (
int eN = 0; eN < nElements_global; eN++)
2998 for (
int j = 0; j < nDOF_trial_element; j++) {
2999 int eN_nDOF_trial_element = eN * nDOF_trial_element;
3000 u_free_dof[r_l2g.data()[eN_nDOF_trial_element + j]] = u_dof.data()[u_l2g.data()[eN_nDOF_trial_element + j]];
3001 u_free_dof_old[r_l2g.data()[eN_nDOF_trial_element + j]] = u_dof_old.data()[u_l2g.data()[eN_nDOF_trial_element + j]];
3003 for (
int i = 0; i < NNZ; i++) {
3004 TransportMatrix[i] = 0.;
3005 TransportMatrixConsistent[i] = 0.;
3006 TransportMatrixn[i] = 0.;
3007 TransportMatrixConsistentn[i] = 0.;
3012 for (
int eN = 0; eN < nElements_global; eN++) {
3013 const int eN_nDOF_trial_element = eN * nDOF_trial_element;
3014 for (
int k = 0; k < nQuadraturePoints_element; k++) {
3015 const int eN_k = eN * nQuadraturePoints_element + k;
3016 double jac[nSpace * nSpace], jacDet, jacInv[nSpace * nSpace], x, y,
z;
3017 ck.calculateMapping_element(eN, k, mesh_dof.data(), mesh_l2g.data(),
3018 mesh_trial_ref.data(), mesh_grad_trial_ref.data(),
3019 jac, jacDet, jacInv, x, y,
z);
3020 const double dV = fabs(jacDet) * dV_ref.data()[k];
3021 for (
int i = 0; i < nDOF_test_element; i++) {
3022 const int eN_i = eN * nDOF_test_element + i;
3023 const int free_gi = r_l2g.data()[eN_i];
3024 const double u_test_dV = u_test_ref.data()[k * nDOF_trial_element + i] * dV;
3025 rho_dof[free_gi] += q_rho.data()[eN_k] * u_test_dV;
3026 ML_rho[free_gi] += u_test_dV;
3030 for (
int i = 0; i < numDOFs; ++i) {
3031 if (ML_rho[i] > 0.0) rho_dof[i] /= ML_rho[i];
3032 else rho_dof[i] = rho;
3042 std::vector<double> rho_n_phi_dof(numDOFs_n, 0.0);
3044 std::vector<double> rho_n_phi_dof_old(numDOFs_n, 0.0);
3047 std::vector<double> dphiN_dz_dof(numDOFs_n, 0.0);
3048 std::vector<double> dphiN_dp_dof(numDOFs_n, 0.0);
3049 std::vector<double> rho_w_dof(numDOFs_n, 0.0);
3050 std::vector<double> rho_n_dof(numDOFs_n, 0.0);
3051 std::vector<double> pc_dof(numDOFs_n, 0.0);
3052 std::vector<double> dpc_dof(numDOFs_n, 0.0);
3053 std::vector<double> krn_dof(numDOFs_n, 0.0);
3054 std::vector<double> dkrn_dof(numDOFs_n, 0.0);
3055 std::vector<double> rho_n_dof_old(numDOFs_n, 0.0);
3056 std::vector<double> pc_dof_old(numDOFs_n, 0.0);
3057 std::vector<double> krn_dof_old(numDOFs_n, 0.0);
3060 std::vector<double> pc_uncap_dof(numDOFs_n, 0.0);
3061 std::vector<double> dpc_uncap_dof(numDOFs_n, 0.0);
3062 std::vector<double> pc_uncap_dof_old(numDOFs_n, 0.0);
3063 std::vector<double> ML_n(numDOFs_n, 0.0);
3064 std::vector<double> Sg_dof_old(numDOFs_n, 0.0);
3065 std::vector<double> X_dof_old(numDOFs_n, 0.0);
3066 for (
int eN = 0; eN < nElements_global; eN++) {
3067 const int mat_eN_proj = elementMaterialTypes.data()[eN];
3068 const double phi_eN = thetaR.data()[mat_eN_proj] + thetaSR.data()[mat_eN_proj];
3069 const double alpha_eN_p = alpha.data()[mat_eN_proj];
3070 const double n_vg_eN_p =
n.data()[mat_eN_proj];
3071 const double krn_end_p = krn_end.data()[mat_eN_proj];
3072 const double S_wr_p = thetaR.data()[mat_eN_proj] / phi_eN;
3073 const double one_m_Sr_p = 1.0 - S_wr_p;
3074 const double Se_trap_L3043 = 1.0 - S_gr.data()[mat_eN_proj] / one_m_Sr_p;
3075 const int eN_nDOF_trial_element = eN * nDOF_trial_element;
3076 for (
int k = 0; k < nQuadraturePoints_element; k++) {
3077 double jac[nSpace * nSpace], jacDet, jacInv[nSpace * nSpace], x_p, y_p, z_p;
3078 ck.calculateMapping_element(eN, k, mesh_dof.data(), mesh_l2g.data(),
3079 mesh_trial_ref.data(), mesh_grad_trial_ref.data(),
3080 jac, jacDet, jacInv, x_p, y_p, z_p);
3081 const double dV = std::fabs(jacDet) * dV_ref.data()[k];
3083 double u_w_p = 0.0, u_n_p = 0.0;
3084 ck.valFromDOF(u_dof.data(),
3085 &u_l2g.data()[eN_nDOF_trial_element],
3086 &u_trial_ref.data()[k * nDOF_trial_element], u_w_p);
3087 ck.valFromDOF(u_dof_n.data(),
3088 &u_l2g_n.data()[eN_nDOF_trial_element],
3089 &u_trial_ref.data()[k * nDOF_trial_element], u_n_p);
3091 const double Se_p_raw = (1.0 - u_n_p - S_wr_p) / one_m_Sr_p;
3092 double Se_p, dSe_du_n_p;
3093 if (Se_p_raw <= 0.0) { Se_p = 0.0; dSe_du_n_p = 0.0; }
3094 else if (Se_p_raw >= 1.0) { Se_p = 1.0; dSe_du_n_p = 0.0; }
3095 else { Se_p = Se_p_raw; dSe_du_n_p = -1.0 / one_m_Sr_p; }
3096 double pc_p = 0.0, dpc_dSe_p = 0.0, d2pc_p_unused = 0.0;
3101 const double dpc_dSn_p = dpc_dSe_p * dSe_du_n_p;
3102 double krn_p = 0.0, dkrn_dSe_p = 0.0;
3107 krn_p = krn_p * krn_end_p / mu_n;
3108 const double dkrn_dSn_p = dkrn_dSe_p * krn_end_p * dSe_du_n_p / mu_n;
3109 const double rho_n_p = rho_n_compressible ? (rho_n * exp(fmin((u_w_p + pc_p) * inv_p_ref_n, 50.0))) : rho_n;
3110 const double phi_rho_n_qp = phi_eN * rho_n_p;
3112 double u_w_p_old = 0.0, u_n_p_old = 0.0;
3113 ck.valFromDOF(u_dof_old.data(),
3114 &u_l2g.data()[eN_nDOF_trial_element],
3115 &u_trial_ref.data()[k * nDOF_trial_element], u_w_p_old);
3116 ck.valFromDOF(u_dof_n_old.data(),
3117 &u_l2g_n.data()[eN_nDOF_trial_element],
3118 &u_trial_ref.data()[k * nDOF_trial_element], u_n_p_old);
3119 const double Se_p_old_raw = (1.0 - u_n_p_old - S_wr_p) / one_m_Sr_p;
3121 if (Se_p_old_raw <= 0.0) Se_p_old = 0.0;
3122 else if (Se_p_old_raw >= 1.0) Se_p_old = 1.0;
3123 else Se_p_old = Se_p_old_raw;
3124 double pc_p_old = 0.0, dpc_p_old_unused = 0.0, d2pc_p_old_unused = 0.0;
3129 double krn_p_old = 0.0, dkrn_p_old_unused = 0.0;
3134 krn_p_old = krn_p_old * krn_end_p / mu_n;
3135 const double rho_n_p_old = rho_n_compressible ? (rho_n * exp(fmin((u_w_p_old + pc_p_old) * inv_p_ref_n, 50.0))) : rho_n;
3137 const int eN_k_proj = eN * nQuadraturePoints_element + k;
3138 const double rho_w_qp_proj = q_rho.data()[eN_k_proj];
3139 const double z_cl_pr = fmin(fmax(u_n_p, 1.0e-8), 1.0 - 1.0e-8);
3140 const double p_cl_pr = fmax(u_w_p, 1.0e2);
3141 const double z_cl_pr_old = fmin(fmax(u_n_p_old, 1.0e-8), 1.0 - 1.0e-8);
3142 const double p_cl_pr_old = fmax(u_w_p_old, 1.0e2);
3147 const double Sa_pr = 1.0 - fs_pr.
S_g;
3148 const double N_pr = fs_pr.
rho_g*fs_pr.
S_g + fs_pr.
rho_a*Sa_pr;
3153 const double N_pr_old = fs_pr_old.
rho_g*fs_pr_old.
S_g
3154 + fs_pr_old.
rho_a*(1.0 - fs_pr_old.
S_g);
3155 const double phiN_qp = phi_eN * N_pr;
3156 const double phiN_old_qp = phi_eN * N_pr_old;
3157 const double dphiN_dz_qp = phi_eN * dN_dz_pr;
3158 const double dphiN_dp_qp = phi_eN * dN_dp_pr;
3159 for (
int i = 0; i < nDOF_test_element; i++) {
3160 const int eN_i = eN * nDOF_test_element + i;
3161 const int gi = u_l2g_n.data()[eN_i];
3162 const double u_test_dV = u_test_ref.data()[k * nDOF_trial_element + i] * dV;
3164 rho_n_phi_dof[gi] += phiN_qp * u_test_dV;
3165 rho_n_phi_dof_old[gi] += phiN_old_qp * u_test_dV;
3166 dphiN_dz_dof[gi] += dphiN_dz_qp * u_test_dV;
3167 dphiN_dp_dof[gi] += dphiN_dp_qp * u_test_dV;
3169 rho_w_dof[gi] += rho_w_qp_proj * u_test_dV;
3170 rho_n_dof[gi] += rho_n_p * u_test_dV;
3171 pc_dof[gi] += pc_p * u_test_dV;
3172 dpc_dof[gi] += dpc_dSn_p * u_test_dV;
3173 krn_dof[gi] += krn_p * u_test_dV;
3174 dkrn_dof[gi] += dkrn_dSn_p * u_test_dV;
3175 rho_n_dof_old[gi] += rho_n_p_old * u_test_dV;
3176 pc_dof_old[gi] += pc_p_old * u_test_dV;
3177 krn_dof_old[gi] += krn_p_old * u_test_dV;
3179 Sg_dof_old[gi] += fs_pr_old.
S_g * u_test_dV;
3180 X_dof_old[gi] += fs_pr_old.
X * u_test_dV;
3181 ML_n[gi] += u_test_dV;
3188 for (
int i = 0; i < numDOFs_n; ++i) {
3189 if (ML_n[i] > 0.0) {
3190 rho_n_phi_dof[i] /= ML_n[i];
3191 rho_n_phi_dof_old[i] /= ML_n[i];
3192 dphiN_dz_dof[i] /= ML_n[i];
3193 dphiN_dp_dof[i] /= ML_n[i];
3194 rho_w_dof[i] /= ML_n[i];
3195 rho_n_dof[i] /= ML_n[i];
3196 pc_dof[i] /= ML_n[i];
3197 dpc_dof[i] /= ML_n[i];
3198 krn_dof[i] /= ML_n[i];
3199 dkrn_dof[i] /= ML_n[i];
3200 rho_n_dof_old[i] /= ML_n[i];
3201 pc_dof_old[i] /= ML_n[i];
3202 krn_dof_old[i] /= ML_n[i];
3203 Sg_dof_old[i] /= ML_n[i];
3204 X_dof_old[i] /= ML_n[i];
3206 rho_n_phi_dof[i] = thetaR.data()[0] + thetaSR.data()[0];
3207 rho_n_phi_dof_old[i] = thetaR.data()[0] + thetaSR.data()[0];
3209 rho_n_dof[i] = rho_n;
3210 rho_n_dof_old[i] = rho_n;
3214 rho_n_phi_dof[i] = std::max(rho_n_phi_dof[i], 1.0e-16);
3219 double psi[numDOFs], eta[numDOFs], global_entropy_residual[numDOFs], boundary_integral[numDOFs];
3220 for (
int i = 0; i < numDOFs; i++) {
3224 double solni = 1.0 * u_free_dof_old[i];
3226 global_entropy_residual[i] = 0.;
3228 boundary_integral[i] = 0.;
3241 for (
int eN = 0; eN < nElements_global; eN++) {
3242 const int eN_nDOF_trial_element = eN * nDOF_trial_element;
3243 const int eN_nDOF_mesh_trial_element = eN * nDOF_mesh_trial_element;
3245 double elementResidual_u[nDOF_test_element], element_entropy_residual[nDOF_test_element], Phi[nDOF_trial_element], Phi_n[nDOF_trial_element];
3246 double elementTransport[nDOF_test_element][nDOF_trial_element], elementTransportConsistent[nDOF_test_element][nDOF_trial_element];
3247 double elementTransportn[nDOF_test_element][nDOF_trial_element], elementTransportConsistentn[nDOF_test_element][nDOF_trial_element];
3248 for (
int j = 0; j < nDOF_trial_element; j++) {
3249 const int u_gj = u_l2g.data()[eN_nDOF_trial_element + j];
3250 const int free_gj = r_l2g.data()[eN_nDOF_trial_element + j];
3251 const int x_gj = mesh_l2g.data()[eN_nDOF_mesh_trial_element + j];
3252 const double rho_node_j = rho_dof[free_gj];
3253 Phi[j] = u_dof.data()[u_gj];
3254 Phi_n[j] = u_dof_old.data()[u_gj];
3255 for (
int I = 0; I < nSpace; I++) {
3257 Phi[j] -= rho_node_j * mesh_dof.data()[x_gj * 3 + I] * gravity[I];
3258 Phi_n[j] -= rho_node_j * mesh_dof.data()[x_gj * 3 + I] * gravity[I];
3261 for (
int i = 0; i < nDOF_test_element; i++) {
3262 elementResidual_u[i] = 0.0;
3263 element_entropy_residual[i] = 0.0;
3264 for (
int j = 0; j < nDOF_trial_element; j++) {
3265 elementTransport[i][j] = 0.0;
3266 elementTransportConsistent[i][j] = 0.0;
3267 elementTransportn[i][j] = 0.0;
3268 elementTransportConsistentn[i][j] = 0.0;
3272 for (
int k = 0; k < nQuadraturePoints_element; k++) {
3274 int eN_k = eN * nQuadraturePoints_element + k, eN_k_nSpace = eN_k * nSpace;
3277 aux_entropy_residual = 0.,
3278 DENTROPY_un, DENTROPY_uni,
3280 u = 0.0, un = 0.0, grad_phi[nSpace], grad_phi_n[nSpace], grad_u_velocity[nSpace], velocity_loc[nSpace], u_test_dV[nDOF_trial_element], u_grad_trial[nDOF_trial_element * nSpace], u_grad_test_dV[nDOF_test_element * nSpace],
3282 jac[nSpace * nSpace], jacDet, jacInv[nSpace * nSpace], dV, x, y,
z,
xt, yt, zt, m, dm,
f[nSpace],
df[nSpace], a[
nnz], da[
nnz], as[
nnz], mn, dmn, fn[nSpace], dfn[nSpace], an[
nnz], dan[
nnz], asn[
nnz];
3284 ck.calculateMapping_element(eN, k, mesh_dof.data(), mesh_l2g.data(), mesh_trial_ref.data(), mesh_grad_trial_ref.data(), jac, jacDet, jacInv, x, y,
z);
3285 ck.calculateMappingVelocity_element(eN, k, mesh_velocity_dof.data(), mesh_l2g.data(), mesh_trial_ref.data(),
xt, yt, zt);
3286 dV = fabs(jacDet) * dV_ref.data()[k];
3288 ck.valFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], &u_trial_ref.data()[k * nDOF_trial_element],
u);
3290 ck.valFromDOF(u_dof_old.data(), &u_l2g.data()[eN_nDOF_trial_element], &u_trial_ref.data()[k * nDOF_trial_element], un);
3292 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k * nDOF_trial_element * nSpace], jacInv, u_grad_trial);
3293 ck.gradFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], u_grad_trial, grad_u_velocity);
3301 for (
int I = 0; I < nSpace; I++) {
3303 grad_phi_n[I] = 0.0;
3305 for (
int j = 0; j < nDOF_trial_element; j++) {
3306 u_test_dV[j] = u_test_ref.data()[k * nDOF_trial_element + j] * dV;
3307 for (
int I = 0; I < nSpace; I++) {
3308 grad_phi_n[I] += Phi_n[j] * u_grad_trial[j * nSpace + I];
3309 grad_phi[I] += Phi[j] * u_grad_trial[j * nSpace + I];
3310 u_grad_test_dV[j * nSpace + I] = u_grad_trial[j * nSpace + I] * dV;
3316 double Kr, dKr, Krn, dKrn, thetaW, thetaWn;
3317 const double rho_local = q_rho.data()[eN_k];
3318 const double rho_velocity = std::fabs(rho_local) > 1.0e-12 ? rho_local : rho;
3322 double dm_du_n_qp_n = 0.0, dkr_du_n_qp_n = 0.0;
3323 double dm_du_n_qp = 0.0, dkr_du_n_qp = 0.0;
3324 double df_du_n_qp_n[nSpace], df_du_n_qp[nSpace];
3325 double da_du_n_qp_n[
nnz], da_du_n_qp[
nnz];
3326 for (
int I = 0; I < nSpace; I++) { df_du_n_qp_n[I] = 0.0; df_du_n_qp[I] = 0.0; }
3327 for (
int ii = 0; ii <
nnz; ii++) { da_du_n_qp_n[ii] = 0.0; da_du_n_qp[ii] = 0.0; }
3328 double u_n_qp = 0.0, u_n_qp_old = 0.0;
3329 ck.valFromDOF(u_dof_n.data(),
3330 &u_l2g.data()[eN_nDOF_trial_element],
3331 &u_trial_ref.data()[k * nDOF_trial_element], u_n_qp);
3332 ck.valFromDOF(u_dof_n_old.data(),
3333 &u_l2g.data()[eN_nDOF_trial_element],
3334 &u_trial_ref.data()[k * nDOF_trial_element], u_n_qp_old);
3336 alpha.data()[elementMaterialTypes[eN]],
n.data()[elementMaterialTypes[eN]],
3337 thetaR.data()[elementMaterialTypes[eN]], thetaSR.data()[elementMaterialTypes[eN]],
3338 &KWs.data()[elementMaterialTypes[eN] *
nnz], un, u_n_qp_old,
3339 mn, dmn, dm_du_n_qp_n, fn, dfn, df_du_n_qp_n, an, dan, da_du_n_qp_n,
3340 asn, Krn, dKrn, dkr_du_n_qp_n, thetaWn);
3342 alpha.data()[elementMaterialTypes[eN]],
n.data()[elementMaterialTypes[eN]],
3343 thetaR.data()[elementMaterialTypes[eN]], thetaSR.data()[elementMaterialTypes[eN]],
3344 &KWs.data()[elementMaterialTypes[eN] *
nnz],
u, u_n_qp,
3345 m, dm, dm_du_n_qp,
f,
df, df_du_n_qp, a, da, da_du_n_qp,
3346 as, Kr, dKr, dkr_du_n_qp, thetaW);
3347 q_theta.data()[eN_k] = thetaW;
3351 for (
int I = 0; I < nSpace; ++I) {
3352 q_velocity.data()[eN_k_nSpace + I] = grad_u_velocity[I];
3355 double pressure_gradient[nSpace];
3356 for (
int J = 0; J < nSpace; ++J)
3357 pressure_gradient[J] = grad_u_velocity[J] - rho_velocity * gravity.data()[J];
3359 for (
int I = 0; I < nSpace; ++I) {
3361 for (
int ii = a_rowptr.data()[I]; ii < a_rowptr.data()[I+1]; ++ii) {
3362 const int J = a_colind.data()[ii];
3363 acc += (a[ii] / rho_velocity) * pressure_gradient[J];
3365 velocity.data()[eN_k_nSpace + I] = -acc;
3366 velocity_couple.data()[eN_k_nSpace + I] = -acc;
3407 double mesh_velocity[3];
3408 mesh_velocity[0] =
xt;
3409 mesh_velocity[1] = yt;
3410 mesh_velocity[2] = zt;
3412 for (
int I = 0; I < nSpace; I++) {
3413 f[I] -= MOVING_DOMAIN * m * mesh_velocity[I];
3414 velocity_loc[I] =
df[I] * (2.0 * dm * dm / (dm * dm + fmax(1.0e-16, dm * dm)));
3419 calculateCFL(elementDiameter.data()[eN] / degree_polynomial, velocity_loc, cfl.data()[eN_k]);
3425 for (
int I = 0; I < nSpace; I++) aux_entropy_residual += velocity_loc[I] * grad_phi_n[I];
3431 for (
int i = 0; i < nDOF_test_element; i++) {
3433 int eN_i = eN * nDOF_test_element + i;
3434 ML2[r_l2g.data()[eN_i]] += u_test_dV[i];
3437 double uni = u_dof_old.data()[u_l2g.data()[eN_i]];
3439 element_entropy_residual[i] += (DENTROPY_un - DENTROPY_uni) * aux_entropy_residual * u_test_dV[i];
3442 elementResidual_u[i] += m * u_test_dV[i];
3447 for (
int j = 0; j < nDOF_trial_element; j++) {
3448 int j_nSpace = j * nSpace;
3449 int i_nSpace = i * nSpace;
3450 elementTransport[i][j] +=
ck.SimpleDiffusionJacobian_weak(a_rowptr.data(), a_colind.data(), as, &u_grad_trial[j_nSpace], &u_grad_test_dV[i_nSpace]);
3451 elementTransportConsistent[i][j] +=
ck.SimpleDiffusionJacobian_weak(a_rowptr.data(), a_colind.data(), a, &u_grad_trial[j_nSpace], &u_grad_test_dV[i_nSpace]);
3452 elementTransportn[i][j] +=
ck.SimpleDiffusionJacobian_weak(a_rowptr.data(), a_colind.data(), asn, &u_grad_trial[j_nSpace], &u_grad_test_dV[i_nSpace]);
3453 elementTransportConsistentn[i][j] +=
ck.SimpleDiffusionJacobian_weak(a_rowptr.data(), a_colind.data(), an, &u_grad_trial[j_nSpace], &u_grad_test_dV[i_nSpace]);
3457 q_u.data()[eN_k] =
u;
3458 q_m.data()[eN_k] = m;
3463 for (
int i = 0; i < nDOF_test_element; i++) {
3464 int eN_i = eN * nDOF_test_element + i;
3465 int gi = r_l2g.data()[eN_i];
3468 global_entropy_residual[gi] += element_entropy_residual[i];
3470 for (
int j = 0; j < nDOF_trial_element; j++) {
3471 int eN_i_j = eN_i * nDOF_trial_element + j;
3472 TransportMatrix[csrRowIndeces_CellLoops.data()[eN_i] + csrColumnOffsets_CellLoops.data()[eN_i_j]] += elementTransport[i][j];
3473 TransportMatrixConsistent[csrRowIndeces_CellLoops.data()[eN_i] + csrColumnOffsets_CellLoops.data()[eN_i_j]] += elementTransportConsistent[i][j];
3474 TransportMatrixn[csrRowIndeces_CellLoops.data()[eN_i] + csrColumnOffsets_CellLoops.data()[eN_i_j]] += elementTransportn[i][j];
3475 TransportMatrixConsistentn[csrRowIndeces_CellLoops.data()[eN_i] + csrColumnOffsets_CellLoops.data()[eN_i_j]] += elementTransportConsistentn[i][j];
3504 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++) {
3505 int ebN = exteriorElementBoundariesArray.data()[ebNE], eN = elementBoundaryElementsArray.data()[ebN * 2 + 0], ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN * 2 + 0], eN_nDOF_trial_element = eN * nDOF_trial_element;
3506 double elementResidual_u[nDOF_test_element];
3507 for (
int i = 0; i < nDOF_test_element; i++) { elementResidual_u[i] = 0.0; }
3508 for (
int kb = 0; kb < nQuadraturePoints_elementBoundary; kb++) {
3509 int ebNE_kb = ebNE * nQuadraturePoints_elementBoundary + kb, ebNE_kb_nSpace = ebNE_kb * nSpace, ebN_local_kb = ebN_local * nQuadraturePoints_elementBoundary + kb, ebN_local_kb_nSpace = ebN_local_kb * nSpace;
3510 double u_ext = 0.0, un_ext, grad_u_ext[nSpace], m_ext = 0.0, dm_ext = 0.0, f_ext[nSpace], df_ext[nSpace], a_ext[
nnz], da_ext[
nnz], as_ext[
nnz],
3511 mn_ext = 0.0, dmn_ext = 0.0, fn_ext[nSpace], dfn_ext[nSpace], an_ext[
nnz], dan_ext[
nnz], asn_ext[
nnz], flux_ext = 0.0, bflux_ext = 0.0,
3513 bc_u_ext = 0.0, bc_grad_u_ext[nSpace], bc_m_ext = 0.0, bc_dm_ext = 0.0, bc_f_ext[nSpace], bc_df_ext[nSpace], bc_a_ext[
nnz], bc_da_ext[
nnz], bc_as_ext[
nnz], jac_ext[nSpace * nSpace], jacDet_ext, jacInv_ext[nSpace * nSpace], boundaryJac[nSpace * (nSpace - 1)], metricTensor[(nSpace - 1) * (nSpace - 1)], metricTensorDetSqrt, dS, u_test_dS[nDOF_test_element], u_grad_trial_trace[nDOF_trial_element * nSpace], normal[3], x_ext, y_ext, z_ext, xt_ext, yt_ext, zt_ext, integralScaling, G[nSpace * nSpace], G_dd_G, tr_G, fluxJacobian_u_u[nDOF_trial_element], bfluxJacobian_u_u[nDOF_trial_element], fluxJacobian_un_un[nDOF_trial_element];
3518 ck.calculateMapping_elementBoundary(eN, ebN_local, kb, ebN_local_kb, mesh_dof.data(), mesh_l2g.data(), mesh_trial_trace_ref.data(), mesh_grad_trial_trace_ref.data(), boundaryJac_ref.data(), jac_ext, jacDet_ext, jacInv_ext, boundaryJac, metricTensor, metricTensorDetSqrt,
3519 normal_ref.data(), normal, x_ext, y_ext, z_ext);
3520 ck.calculateMappingVelocity_elementBoundary(eN, ebN_local, kb, ebN_local_kb, mesh_velocity_dof.data(), mesh_l2g.data(), mesh_trial_trace_ref.data(), xt_ext, yt_ext, zt_ext, normal, boundaryJac, metricTensor, integralScaling);
3521 dS = ((1.0 - MOVING_DOMAIN) * metricTensorDetSqrt + MOVING_DOMAIN * integralScaling) * dS_ref.data()[kb];
3524 ck.calculateG(jacInv_ext, G, G_dd_G, tr_G);
3527 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace * nDOF_trial_element], jacInv_ext, u_grad_trial_trace);
3529 ck.valFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], &u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element], u_ext);
3530 ck.valFromDOF(u_dof_old.data(), &u_l2g.data()[eN_nDOF_trial_element], &u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element], un_ext);
3531 ck.gradFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], u_grad_trial_trace, grad_u_ext);
3541 for (
int j = 0; j < nDOF_trial_element; j++) { u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb * nDOF_test_element + j] * dS; }
3545 bc_u_ext = isDOFBoundary_u.data()[ebNE_kb] * ebqe_bc_u_ext.data()[ebNE_kb] + (1 - isDOFBoundary_u.data()[ebNE_kb]) * u_ext;
3549 double bc_Kr, bc_dKr,bc_Kr_ext, bc_dKr_ext, bc_Krn, bc_dKrn, thetaW_ext, thetaWn_ext, thetaW_bc_ext;
3550 const double rho_ext = ebqe_rho.data()[ebNE_kb];
3551 const double rho_velocity_ext = std::fabs(rho_ext) > 1.0e-12 ? rho_ext : rho;
3554 double dm_du_n_ext = 0.0, dkr_du_n_ext = 0.0;
3555 double dmn_du_n_ext = 0.0, dkrn_du_n_ext = 0.0;
3556 double bc_dm_du_n = 0.0, bc_dkr_du_n = 0.0;
3557 double df_du_n_ext[nSpace], dfn_du_n_ext[nSpace], bc_df_du_n[nSpace];
3558 double da_du_n_ext[
nnz], dan_du_n_ext[
nnz], bc_da_du_n[
nnz];
3559 for (
int I = 0; I < nSpace; I++) {
3560 df_du_n_ext[I] = 0.0; dfn_du_n_ext[I] = 0.0; bc_df_du_n[I] = 0.0;
3562 for (
int ii = 0; ii <
nnz; ii++) {
3563 da_du_n_ext[ii] = 0.0; dan_du_n_ext[ii] = 0.0; bc_da_du_n[ii] = 0.0;
3565 double u_n_ext_qp = 0.0, u_n_ext_qp_old = 0.0;
3566 ck.valFromDOF(u_dof_n.data(), &u_l2g.data()[eN_nDOF_trial_element],
3567 &u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element], u_n_ext_qp);
3568 ck.valFromDOF(u_dof_n_old.data(), &u_l2g.data()[eN_nDOF_trial_element],
3569 &u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element], u_n_ext_qp_old);
3570 const double bc_u_n_ext_qp = isDOFBoundary_n.data()[ebNE_kb] * ebqe_bc_u_n_ext.data()[ebNE_kb]
3571 + (1 - isDOFBoundary_n.data()[ebNE_kb]) * u_n_ext_qp;
3573 alpha.data()[elementMaterialTypes.data()[eN]],
n.data()[elementMaterialTypes.data()[eN]],
3574 thetaR.data()[elementMaterialTypes.data()[eN]], thetaSR.data()[elementMaterialTypes.data()[eN]],
3575 &KWs.data()[elementMaterialTypes.data()[eN] *
nnz], u_ext, u_n_ext_qp,
3576 m_ext, dm_ext, dm_du_n_ext, f_ext, df_ext, df_du_n_ext, a_ext, da_ext, da_du_n_ext,
3577 as_ext, bc_Kr, bc_dKr, dkr_du_n_ext, thetaW_ext);
3579 alpha.data()[elementMaterialTypes.data()[eN]],
n.data()[elementMaterialTypes.data()[eN]],
3580 thetaR.data()[elementMaterialTypes.data()[eN]], thetaSR.data()[elementMaterialTypes.data()[eN]],
3581 &KWs.data()[elementMaterialTypes.data()[eN] *
nnz], un_ext, u_n_ext_qp_old,
3582 mn_ext, dmn_ext, dmn_du_n_ext, fn_ext, dfn_ext, dfn_du_n_ext, an_ext, dan_ext, dan_du_n_ext,
3583 asn_ext, bc_Krn, bc_dKrn, dkrn_du_n_ext, thetaWn_ext);
3585 alpha.data()[elementMaterialTypes.data()[eN]],
n.data()[elementMaterialTypes.data()[eN]],
3586 thetaR.data()[elementMaterialTypes.data()[eN]], thetaSR.data()[elementMaterialTypes.data()[eN]],
3587 &KWs.data()[elementMaterialTypes.data()[eN] *
nnz], bc_u_ext, bc_u_n_ext_qp,
3588 bc_m_ext, bc_dm_ext, bc_dm_du_n, bc_f_ext, bc_df_ext, bc_df_du_n, bc_a_ext, bc_da_ext, bc_da_du_n,
3589 bc_as_ext, bc_Kr_ext, bc_dKr_ext, bc_dkr_du_n, thetaW_bc_ext);
3590 ebqe_theta.data()[ebNE_kb] = thetaW_ext;
3604 double ext_pressure_gradient[nSpace];
3605 for (
int J = 0; J < nSpace; ++J)
3606 ext_pressure_gradient[J] = grad_u_ext[J] - rho_velocity_ext * gravity.data()[J];
3608 for (
int I = 0; I < nSpace; ++I) {
3610 for (
int ii = a_rowptr.data()[I]; ii < a_rowptr.data()[I+1]; ++ii) {
3611 const int J = a_colind.data()[ii];
3612 acc += (a_ext[ii] / rho_velocity_ext) * ext_pressure_gradient[J];
3614 ebqe_velocity_ext.data()[ebNE_kb_nSpace + I] = -acc;
3615 ebqe_velocity_ext_couple.data()[ebNE_kb_nSpace + I] = -acc;
3617 bool useConsistentFlux=
false;
3618 double grad_u_n_ext_b[nSpace];
3619 ck.gradFromDOF(u_dof_n.data(), &u_l2g.data()[eN_nDOF_trial_element],
3620 u_grad_trial_trace, grad_u_n_ext_b);
3621 const int mat_b0 = elementMaterialTypes.data()[eN];
3622 const double alpha_b0 = alpha.data()[mat_b0];
3623 const double n_vg_b0 =
n.data()[mat_b0];
3624 const double krn_end_b0 = krn_end.data()[mat_b0];
3625 const double *KWs_b0 = &KWs.data()[mat_b0 *
nnz];
3626 const double phi_b0 = thetaR.data()[mat_b0] + thetaSR.data()[mat_b0];
3627 const double S_wr_b0 = thetaR.data()[mat_b0] / phi_b0;
3628 const double one_m_Sr_b0 = 1.0 - S_wr_b0;
3629 const double Se_trap_L3616 = 1.0 - S_gr.data()[mat_b0] / one_m_Sr_b0;
3630 const double z_clb0 = fmin(fmax(u_n_ext_qp, 1.0e-8), 1.0 - 1.0e-8);
3631 const double p_clb0 = fmax(u_ext, 1.0e2);
3634 const double Sab0 = 1.0 - fsb0.
S_g;
3635 const double Se_rawb0 = (Sab0 - S_wr_b0)/one_m_Sr_b0;
3636 double Se_ab0, dSeb0_dp, dSeb0_dz;
3637 if (Se_rawb0<=0.0){Se_ab0=0.0;dSeb0_dp=0.0;dSeb0_dz=0.0;}
3638 else if (Se_rawb0>=1.0){Se_ab0=1.0;dSeb0_dp=0.0;dSeb0_dz=0.0;}
3639 else {Se_ab0=Se_rawb0;dSeb0_dp=-fsb0.
dS_g_dp/one_m_Sr_b0;dSeb0_dz=-fsb0.
dS_g_dz/one_m_Sr_b0;}
3640 double KWrb0=0,DKWrb0=0,thWb0=0,DthWb0=0,KNrb0=0,DKNrb0=0,pcb0=0,dpc_dSeb0=0,d2pcb0=0;
3650 KNrb0 *= krn_end_b0; DKNrb0 *= krn_end_b0;
3651 const double pcpb0 = dpc_dSeb0/one_m_Sr_b0;
3652 const double dpcpb0_dp = (d2pcb0/one_m_Sr_b0)*dSeb0_dp;
3653 const double dpcpb0_dz = (d2pcb0/one_m_Sr_b0)*dSeb0_dz;
3657 const double rgmb0 = fsb0.
rho_g*Mbar_gb0, ramb0 = fsb0.
rho_a*Mbar_ab0;
3659 const double drgmb0_dz = fsb0.
rho_g*fsb0.
dY_dz*dMmb0;
3662 const double Agb0 = fsb0.
rho_g*(1.0-fsb0.
Y), Aab0 = fsb0.
rho_a*(1.0-fsb0.
X);
3664 const double dAgb0_dz = - fsb0.
rho_g*fsb0.
dY_dz;
3668 double valp_b[nSpace], valz_b[nSpace];
3669 for (
int I=0;I<nSpace;I++){
3670 double ugI=0,uaI=0,dugp=0,dugz=0,duap=0,duaz=0;
3671 for (
int ii=a_rowptr.data()[I];ii<a_rowptr.data()[I+1];ii++){
3672 const int J=a_colind.data()[ii];
3673 const double Kii=KWs_b0[ii];
3674 const double Mob_g=KNrb0*Kii/mu_n, Mob_a=KWrb0*Kii;
3675 const double dMobg_dp=(DKNrb0*Kii/mu_n)*dSeb0_dp, dMobg_dz=(DKNrb0*Kii/mu_n)*dSeb0_dz;
3676 const double dMoba_dp=(DKWrb0*Kii)*dSeb0_dp, dMoba_dz=(DKWrb0*Kii)*dSeb0_dz;
3677 const double gJ=gravity.data()[J];
3678 const double gradSa=-(fsb0.
dS_g_dp*grad_u_ext[J]+fsb0.
dS_g_dz*grad_u_n_ext_b[J]);
3679 const double gp_a=grad_u_ext[J]-ramb0*gJ;
3680 const double gp_g=grad_u_ext[J]+pcpb0*gradSa-rgmb0*gJ;
3681 ugI-=Mob_g*gp_g; uaI-=Mob_a*gp_a;
3682 const double dgradSa_dp=-(fsb0.
d2S_g_dp2*grad_u_ext[J]+fsb0.
d2S_g_dpdz*grad_u_n_ext_b[J]);
3683 const double dgradSa_dz=-(fsb0.
d2S_g_dpdz*grad_u_ext[J]+fsb0.
d2S_g_dz2*grad_u_n_ext_b[J]);
3684 const double dgpg_dp=dpcpb0_dp*gradSa+pcpb0*dgradSa_dp-drgmb0_dp*gJ;
3685 const double dgpg_dz=dpcpb0_dz*gradSa+pcpb0*dgradSa_dz-drgmb0_dz*gJ;
3686 dugp-=dMobg_dp*gp_g+Mob_g*dgpg_dp;
3687 dugz-=dMobg_dz*gp_g+Mob_g*dgpg_dz;
3688 duap-=dMoba_dp*gp_a+Mob_a*(-dramb0_dp*gJ);
3689 duaz-=dMoba_dz*gp_a+Mob_a*(-dramb0_dz*gJ);
3691 F0n_b += (Agb0*ugI + Aab0*uaI) * normal[I];
3692 valp_b[I] = (dAgb0_dp*ugI + Agb0*dugp + dAab0_dp*uaI + Aab0*duap);
3693 valz_b[I] = (dAgb0_dz*ugI + Agb0*dugz + dAab0_dz*uaI + Aab0*duaz);
3695 const double penb0 = ebqe_penalty_ext.data()[ebNE_kb];
3698 const int isSeep = isSeepageFace.data()[ebNE];
3699 if (isSeep || isDOFBoundary_u.data()[ebNE_kb]) {
3700 const double bc_u_pen = isSeep ? 0.0 : bc_u_ext;
3701 const double pen_term = penb0*(u_ext - bc_u_pen);
3702 flux_ext = F0n_b + pen_term;
3703 bflux_ext = pen_term;
3704 if (isSeep && flux_ext <= 0.0) { flux_ext = 0.0; bflux_ext = 0.0; }
3706 flux_ext = ebqe_bc_flux_ext[ebNE_kb];
3707 bflux_ext = ebqe_bc_flux_ext[ebNE_kb];
3711 ebqe_flux.data()[ebNE_kb] = flux_ext;
3713 anb_seepage_flux =
seepagefluxcalculator(anb_seepage_flux, isSeepageFace.data()[ebNE], dS, flux_ext);
3714 anb_seepage_flux_n.data()[0] = anb_seepage_flux;
3715 ebqe_u.data()[ebNE_kb] = u_ext;
3719 for (
int i = 0; i < nDOF_test_element; i++) {
3720 if (useConsistentFlux) {
3721 elementResidual_u[i] +=
ck.ExteriorElementBoundaryFlux(flux_ext, u_test_dS[i]);
3723 elementResidual_u[i] +=
ck.ExteriorElementBoundaryFlux(bflux_ext, u_test_dS[i]);
3726 for (
int j = 0; j < nDOF_trial_element; j++) {
3727 if (useConsistentFlux) {
3728 exteriorNumericalFluxJacobian(a_rowptr.data(), a_colind.data(), isDOFBoundary_u.data()[ebNE_kb], normal, a_ext, da_ext, grad_u_ext, &u_grad_trial_trace[j * nSpace], df_ext, u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element + j],
3729 ebqe_penalty_ext.data()[ebNE_kb],
3730 fluxJacobian_u_u[j]);
3737 if (isDOFBoundary_u.data()[ebNE_kb]) {
3738 const double trial_j = u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element + j];
3740 for (
int I=0;I<nSpace;I++){
3741 accp += normal[I]*valp_b[I]*trial_j;
3742 for (
int ii=a_rowptr.data()[I];ii<a_rowptr.data()[I+1];ii++){
3743 const int J=a_colind.data()[ii];
3744 const double Kii=KWs_b0[ii];
3745 const double Mob_g=KNrb0*Kii/mu_n, Mob_a=KWrb0*Kii;
3746 const double dFdgp = Agb0*(-Mob_g*(1.0 - pcpb0*fsb0.
dS_g_dp)) + Aab0*(-Mob_a);
3747 accp += normal[I]*dFdgp*u_grad_trial_trace[j*nSpace + J];
3750 fluxJacobian_u_u[j] = accp;
3751 bfluxJacobian_u_u[j] = penb0*trial_j;
3753 fluxJacobian_u_u[j] = 0.0;
3754 bfluxJacobian_u_u[j] = 0.0;
3758 isDOFBoundary_u.data()[ebNE_kb], normal,
3759 asn_ext, dan_ext, grad_u_ext,
3760 &u_grad_trial_trace[j * nSpace], dfn_ext,
3761 u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element + j],
3762 ebqe_penalty_ext.data()[ebNE_kb],
3763 fluxJacobian_un_un[j]);
3768 for (
int i = 0; i < nDOF_test_element; i++) {
3769 int eN_i = eN * nDOF_test_element + i;
3770 for (
int j = 0; j < nDOF_trial_element; j++) {
3772 if (useConsistentFlux) {
3773 globalJacobian.data()[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += fluxJacobian_u_u[j] * u_test_dS[i];
3775 globalJacobian.data()[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += bfluxJacobian_u_u[j] * u_test_dS[i];
3776 TransportMatrix[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += fluxJacobian_u_u[j] * u_test_dS[i];
3777 TransportMatrixConsistent[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += fluxJacobian_u_u[j] * u_test_dS[i];
3778 TransportMatrixn[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += fluxJacobian_un_un[j] * u_test_dS[i];
3779 TransportMatrixConsistentn[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += fluxJacobian_un_un[j] * u_test_dS[i];
3784 for (
int i = 0; i < nDOF_test_element; i++) {
3785 int eN_i = eN * nDOF_test_element + i;
3786 globalResidual.data()[offset_u + stride_u * u_l2g.data()[eN_i]] += elementResidual_u[i];
3793 std::vector<double> cflux(numDOFs, 0.0);
3794 for (
int i = 0; i < numDOFs_u; i++) {
3795 double gi[nSpace], Cij[nSpace], xi[nSpace], etaMaxi, etaMini;
3796 const int node_i = freeDOFToNode_u.data()[i];
3797 double solni = u_free_dof_old[i];
3798 for (
int I = 0; I < nSpace; I++) {
3799 solni -= rho_dof[i] * gravity.data()[I] * mesh_dof.data()[node_i * 3 + I];
3804 etaMaxi = fabs(eta[i]);
3805 etaMini = fabs(eta[i]);
3808 for (
int I = 0; I < nSpace; I++) {
3810 xi[I] = mesh_dof.data()[node_i * 3 + I];
3813 double alpha_numerator_pos = 0., alpha_numerator_neg = 0., alpha_denominator_pos = 0., alpha_denominator_neg = 0.;
3814 for (
int offset = csrRowIndeces_DofLoops.data()[i]; offset < csrRowIndeces_DofLoops.data()[i + 1]; offset++) {
3815 int j = csrColumnOffsets_DofLoops.data()[offset];
3816 const int full_offset = full_offset_from_compact(i, j);
3817 assert(full_offset >= 0);
3818 const int node_j = freeDOFToNode_u.data()[j];
3822 etaMaxi = fmax(etaMaxi, fabs(eta[j]));
3823 etaMini = fmin(etaMini, fabs(eta[j]));
3825 double solnj = u_free_dof_old[j];
3826 for (
int I = 0; I < nSpace; I++) {
3827 solnj -= rho_dof[j] * gravity.data()[I] * mesh_dof.data()[node_j * 3 + I];
3830 Cij[0] = Cx[full_offset];
3832 Cij[1] = Cy[full_offset];
3835 Cij[2] = Cz[full_offset];
3838 for (
int I = 0; I < nSpace; I++) gi[I] += Cij[I] * solnj;
3841 double alpha_num = solni - solnj;
3842 if (alpha_num >= 0.) {
3843 alpha_numerator_pos += alpha_num;
3844 alpha_denominator_pos += alpha_num;
3846 alpha_numerator_neg += alpha_num;
3847 alpha_denominator_neg += fabs(alpha_num);
3853 for (
int I = 0; I < nSpace; I++) gi[I] /= ML.data()[i];
3857 global_entropy_residual[i] *= etaMini == etaMaxi ? 0. : 2 *
cE / (etaMaxi - etaMini);
3858 quantDOFs.data()[i] = fabs(global_entropy_residual[i]);
3862 double SumPos = 0., SumNeg = 0.;
3863 for (
int offset = csrRowIndeces_DofLoops.data()[i]; offset < csrRowIndeces_DofLoops.data()[i + 1]; offset++) {
3864 int j = csrColumnOffsets_DofLoops.data()[offset];
3865 const int full_offset = full_offset_from_compact(i, j);
3866 assert(full_offset >= 0);
3868 double gi_times_x = 0.;
3869 for (
int I = 0; I < nSpace; I++) {
3870 gi_times_x += gi[I] * delta_x_ij.data()[full_offset * 3 + I];
3873 SumPos += gi_times_x > 0 ? gi_times_x : 0;
3874 SumNeg += gi_times_x < 0 ? gi_times_x : 0;
3876 double sigmaPosi = fmin(1., (fabs(SumNeg) + 1E-15) / (SumPos + 1E-15));
3877 double sigmaNegi = fmin(1., (SumPos + 1E-15) / (fabs(SumNeg) + 1E-15));
3878 double alpha_numi = fabs(sigmaPosi * alpha_numerator_pos + sigmaNegi * alpha_numerator_neg);
3879 double alpha_deni = sigmaPosi * alpha_denominator_pos + sigmaNegi * alpha_denominator_neg;
3881 alpha_numi = fabs(alpha_numerator_pos + alpha_numerator_neg);
3882 alpha_deni = alpha_denominator_pos + alpha_denominator_neg;
3884 double alphai = alpha_numi / (alpha_deni + 1E-15);
3885 quantDOFs.data()[i] = alphai;
3890 std::vector<double> w_lam_g(numDOFs_u, 0.0), w_lam_a(numDOFs_u, 0.0);
3891 std::vector<double> w_dlam_g_dp(numDOFs_u, 0.0), w_dlam_g_dz(numDOFs_u, 0.0);
3892 std::vector<double> w_dlam_a_dp(numDOFs_u, 0.0), w_dlam_a_dz(numDOFs_u, 0.0);
3893 std::vector<double> w_pc(numDOFs_u, 0.0), w_dpc_dp(numDOFs_u, 0.0), w_dpc_dz(numDOFs_u, 0.0);
3894 std::vector<double> w_rgm(numDOFs_u, 0.0), w_ram(numDOFs_u, 0.0);
3895 std::vector<double> w_drgm_dp(numDOFs_u, 0.0), w_drgm_dz(numDOFs_u, 0.0);
3896 std::vector<double> w_dram_dp(numDOFs_u, 0.0), w_dram_dz(numDOFs_u, 0.0);
3897 std::vector<double> w_lam_g_old(numDOFs_u, 0.0), w_lam_a_old(numDOFs_u, 0.0);
3898 std::vector<double> w_pc_old(numDOFs_u, 0.0);
3899 std::vector<double> w_rgm_old(numDOFs_u, 0.0), w_ram_old(numDOFs_u, 0.0);
3902 for (
int i = 0; i < numDOFs_u; i++) {
3903 const int node_i = freeDOFToNode_u.data()[i];
3906 const int z_i = (split_z != 0) ? node2zdof.data()[node_i] : node_i;
3907 const int mat_i = freeDOFMaterialTypes.data()[i];
3908 const double alpha_i = alpha.data()[mat_i];
3909 const double n_vg_i =
n.data()[mat_i];
3910 const double krn_end_i= krn_end.data()[mat_i];
3911 const double phi_i = thetaR.data()[mat_i] + thetaSR.data()[mat_i];
3912 const double S_wr_i = thetaR.data()[mat_i] / phi_i;
3913 const double one_m_Sr_i = 1.0 - S_wr_i;
3914 const double Se_trap_L3923 = 1.0 - S_gr.data()[mat_i] / one_m_Sr_i;
3915 const double cg_i = krn_end_i / mu_n;
3917 const double z_cl = fmin(fmax(u_dof_n.data()[z_i], 1.0e-8), 1.0 - 1.0e-8);
3918 const double p_cl = fmax(u_free_dof[i], 1.0e2);
3921 const double Sa = 1.0 -
f.S_g;
3922 const double Se_raw = (Sa - S_wr_i) / one_m_Sr_i;
3923 double Se, dSe_dp, dSe_dz;
3924 if (Se_raw <= 0.0) { Se = 0.0; dSe_dp = 0.0; dSe_dz = 0.0; }
3925 else if (Se_raw >= 1.0) { Se = 1.0; dSe_dp = 0.0; dSe_dz = 0.0; }
3926 else { Se = Se_raw; dSe_dp = -
f.dS_g_dp/one_m_Sr_i; dSe_dz = -
f.dS_g_dz/one_m_Sr_i; }
3927 double krn=0,dkrn=0,krw=0,dkrw=0,thW=0,DthW=0,pc=0,dpc_dSe=0,d2pc=0;
3937 w_pc[i] = pc; w_dpc_dp[i] = dpc_dSe*dSe_dp; w_dpc_dz[i] = dpc_dSe*dSe_dz;
3940 w_rgm[i] =
f.rho_g*Mg; w_ram[i] =
f.rho_a*Ma;
3941 w_drgm_dp[i] =
f.drho_g_dp*Mg +
f.rho_g*
f.dY_dp*dMm_w; w_drgm_dz[i] =
f.rho_g*
f.dY_dz*dMm_w;
3942 w_dram_dp[i] =
f.drho_a_dp*Ma +
f.rho_a*
f.dX_dp*dMm_w; w_dram_dz[i] =
f.drho_a_dz*Ma +
f.rho_a*
f.dX_dz*dMm_w;
3944 const double yw = 1.0 -
f.Y, dyw_dp = -
f.dY_dp, dyw_dz = -
f.dY_dz;
3945 const double xw = 1.0 -
f.X, dxw_dp = -
f.dX_dp, dxw_dz = -
f.dX_dz;
3946 w_lam_g[i] = cg_i*
f.rho_g*yw*krn;
3947 w_dlam_g_dp[i] = cg_i*(
f.drho_g_dp*yw*krn +
f.rho_g*dyw_dp*krn +
f.rho_g*yw*dkrn*dSe_dp);
3948 w_dlam_g_dz[i] = cg_i*(
f.rho_g*dyw_dz*krn +
f.rho_g*yw*dkrn*dSe_dz);
3949 w_lam_a[i] =
f.rho_a*xw*krw;
3950 w_dlam_a_dp[i] =
f.drho_a_dp*xw*krw +
f.rho_a*dxw_dp*krw +
f.rho_a*xw*dkrw*dSe_dp;
3951 w_dlam_a_dz[i] =
f.drho_a_dz*xw*krw +
f.rho_a*dxw_dz*krw +
f.rho_a*xw*dkrw*dSe_dz;
3953 const double z_cl_o = fmin(fmax(u_dof_n_old.data()[z_i], 1.0e-8), 1.0 - 1.0e-8);
3954 const double p_cl_o = fmax(u_free_dof_old[i], 1.0e2);
3957 const double Sa_o = 1.0 - fo.
S_g;
3958 const double Se_o_raw = (Sa_o - S_wr_i) / one_m_Sr_i;
3959 const double Se_o = Se_o_raw <= 0.0 ? 0.0 : (Se_o_raw >= 1.0 ? 1.0 : Se_o_raw);
3960 double krn_o=0,dkrn_o=0,krw_o=0,dkrw_o=0,thW_o=0,DthW_o=0,pc_o=0,dpc_o=0,d2pc_o=0;
3973 w_rgm_old[i] = fo.
rho_g*Mg_o; w_ram_old[i] = fo.
rho_a*Ma_o;
3974 w_lam_g_old[i] = cg_i*fo.
rho_g*(1.0-fo.
Y)*krn_o;
3975 w_lam_a_old[i] = fo.
rho_a*(1.0-fo.
X)*krw_o;
3984 for (
int i = 0; i < numDOFs_u; i++) {
3986 double sum_abs_dt_times_fH_minus_fL = 0.0, MLi = ML.data()[i];
3987 double Kr, dKr, Krn, dKrn;
3989 double ith_dissipative_term = 0;
3990 double ith_low_order_dissipative_term = 0;
3991 double ith_flux_term = 0;
3992 double ith_consistent_flux_term = 0;
3994 double m, dm,
f[nSpace],
df[nSpace], a[
nnz], da[
nnz], as[
nnz];
3995 double dmn, fn[nSpace], dfn[nSpace], an[
nnz], dan[
nnz], asn[
nnz];
4001 double dm_du_n_fct, dkr_du_n_fct, df_du_n_fct[nSpace], da_du_n_fct[
nnz];
4002 double dmn_du_n_fct, dkrn_du_n_fct, dfn_du_n_fct[nSpace], dan_du_n_fct[
nnz];
4004 const double rho_i = rho_dof[i];
4005 const int node_i = freeDOFToNode_u.data()[i];
4008 const int z_i = (split_z != 0) ? node2zdof.data()[node_i] : node_i;
4010 double thetaW_tmp = 0.0;
4012 for (
int offset = csrRowIndeces_DofLoops.data()[i]; offset < csrRowIndeces_DofLoops.data()[i + 1]; offset++) {
4013 int j = csrColumnOffsets_DofLoops.data()[offset];
4014 const int full_offset = full_offset_from_compact(i, j);
4015 assert(full_offset >= 0);
4016 if (i == j) ii = full_offset;
4017 const double rho_j = rho_dof[j];
4018 const int node_j = freeDOFToNode_u.data()[j];
4019 const double rho_edge = 0.5 * (rho_i + rho_j);
4036 dt_times_fH_minus_fL.data()[full_offset] = 0.0;
4037 if (i == j)
continue;
4038 const double tau = fmax(0.0, -TransportMatrix[full_offset]) / rho_edge;
4039 if (tau == 0.0)
continue;
4040 double g_dot_dx = 0.0;
4041 for (
int I = 0; I < nSpace; I++)
4042 g_dot_dx += gravity.data()[I]
4043 * (mesh_dof.data()[node_j * 3 + I] - mesh_dof.data()[node_i * 3 + I]);
4045 const double rgm_edge = 0.5 * (w_rgm[i] + w_rgm[j]);
4046 const double rgm_edge_old = 0.5 * (w_rgm_old[i] + w_rgm_old[j]);
4047 const double dPhi_g = (u_free_dof[j] + w_pc[j]) - (u_free_dof[i] + w_pc[i])
4048 - rgm_edge * g_dot_dx;
4049 const double dPhi_g_old = (u_free_dof_old[j] + w_pc_old[j]) - (u_free_dof_old[i] + w_pc_old[i])
4050 - rgm_edge_old * g_dot_dx;
4051 const bool up_i_g = (dPhi_g <= 0.0);
4052 const bool up_i_g_old = (dPhi_g_old <= 0.0);
4053 const double lam_g_up = up_i_g ? w_lam_g[i] : w_lam_g[j];
4054 const double lam_g_up_old = up_i_g_old ? w_lam_g_old[i] : w_lam_g_old[j];
4055 const double Fg = Theta * tau * lam_g_up * dPhi_g
4056 + (1.0 - Theta) * tau * lam_g_up_old * dPhi_g_old;
4058 const double ram_edge = 0.5 * (w_ram[i] + w_ram[j]);
4059 const double ram_edge_old = 0.5 * (w_ram_old[i] + w_ram_old[j]);
4060 const double dPhi_a = u_free_dof[j] - u_free_dof[i] - ram_edge * g_dot_dx;
4061 const double dPhi_a_old = u_free_dof_old[j] - u_free_dof_old[i] - ram_edge_old * g_dot_dx;
4062 const bool up_i_a = (dPhi_a <= 0.0);
4063 const bool up_i_a_old = (dPhi_a_old <= 0.0);
4064 const double lam_a_up = up_i_a ? w_lam_a[i] : w_lam_a[j];
4065 const double lam_a_up_old = up_i_a_old ? w_lam_a_old[i] : w_lam_a_old[j];
4066 const double Fa = Theta * tau * lam_a_up * dPhi_a
4067 + (1.0 - Theta) * tau * lam_a_up_old * dPhi_a_old;
4078 const double ddPhig_dpi = -1.0 - w_dpc_dp[i] - 0.5 * w_drgm_dp[i] * g_dot_dx;
4079 const double ddPhig_dzi = - w_dpc_dz[i] - 0.5 * w_drgm_dz[i] * g_dot_dx;
4080 const double ddPhig_dpj = +1.0 + w_dpc_dp[j] - 0.5 * w_drgm_dp[j] * g_dot_dx;
4081 const double ddPhig_dzj = + w_dpc_dz[j] - 0.5 * w_drgm_dz[j] * g_dot_dx;
4082 const double ddPhia_dpi = -1.0 - 0.5 * w_dram_dp[i] * g_dot_dx;
4083 const double ddPhia_dzi = - 0.5 * w_dram_dz[i] * g_dot_dx;
4084 const double ddPhia_dpj = +1.0 - 0.5 * w_dram_dp[j] * g_dot_dx;
4085 const double ddPhia_dzj = - 0.5 * w_dram_dz[j] * g_dot_dx;
4086 const double Tt = Theta * tau;
4087 double dF_dpi = 0.0, dF_dzi = 0.0, dF_dpj = 0.0, dF_dzj = 0.0;
4089 dF_dpi += Tt*lam_g_up*ddPhig_dpi; dF_dzi += Tt*lam_g_up*ddPhig_dzi;
4090 dF_dpj += Tt*lam_g_up*ddPhig_dpj; dF_dzj += Tt*lam_g_up*ddPhig_dzj;
4092 if (up_i_g) { dF_dpi += Tt*w_dlam_g_dp[i]*dPhi_g; dF_dzi += Tt*w_dlam_g_dz[i]*dPhi_g; }
4093 else { dF_dpj += Tt*w_dlam_g_dp[j]*dPhi_g; dF_dzj += Tt*w_dlam_g_dz[j]*dPhi_g; }
4095 dF_dpi += Tt*lam_a_up*ddPhia_dpi; dF_dzi += Tt*lam_a_up*ddPhia_dzi;
4096 dF_dpj += Tt*lam_a_up*ddPhia_dpj; dF_dzj += Tt*lam_a_up*ddPhia_dzj;
4098 if (up_i_a) { dF_dpi += Tt*w_dlam_a_dp[i]*dPhi_a; dF_dzi += Tt*w_dlam_a_dz[i]*dPhi_a; }
4099 else { dF_dpj += Tt*w_dlam_a_dp[j]*dPhi_a; dF_dzj += Tt*w_dlam_a_dz[j]*dPhi_a; }
4103 (void)dF_dpi; (void)dF_dpj; (void)dF_dzi; (void)dF_dzj;
4105 mDotLow.data()[i] = ith_flux_term/MLi;
4106 cflux[i] = ith_consistent_flux_term;
4112 alpha.data()[freeDOFMaterialTypes.data()[i]],
4113 n.data()[freeDOFMaterialTypes.data()[i]], thetaR.data()[freeDOFMaterialTypes.data()[i]], thetaSR.data()[freeDOFMaterialTypes.data()[i]], &KWs.data()[freeDOFMaterialTypes.data()[i] *
nnz],
4114 u_free_dof[i], u_dof_n.data()[z_i],
4115 m, dm, dm_du_n_fct,
f,
df, df_du_n_fct, a, da, da_du_n_fct, as, Kr, dKr, dkr_du_n_fct, thetaW_tmp);
4117 alpha.data()[freeDOFMaterialTypes.data()[i]],
4118 n.data()[freeDOFMaterialTypes.data()[i]], thetaR.data()[freeDOFMaterialTypes.data()[i]], thetaSR.data()[freeDOFMaterialTypes.data()[i]], &KWs.data()[freeDOFMaterialTypes.data()[i] *
nnz],
4119 u_free_dof_old[i], u_dof_n_old.data()[z_i],
4120 mn.data()[i], dmn, dmn_du_n_fct, fn, dfn, dfn_du_n_fct, an, dan, dan_du_n_fct, asn, Krn, dKrn, dkrn_du_n_fct, thetaW_tmp);
4122 globalResidual.data()[offset_u + stride_u * i] += bc_mask.data()[i] * (MLi * (m - mn.data()[i]) / dt - ith_flux_term);
4123 globalJacobian.data()[ii] += bc_mask.data()[i] * (MLi * dm / dt + J_ii) + (1.0 - bc_mask.data()[i]);
4130 const int row_w = offset_u + stride_u * i;
4131 const int col_n = offset_n + stride_n * z_i;
4133 for (
int o = csrRowIndeces_Full.data()[row_w];
4134 o < csrRowIndeces_Full.data()[row_w + 1]; o++) {
4135 if (csrColumnOffsets_Full.data()[o] == col_n) { off_wv = o;
break; }
4138 globalJacobian.data()[off_wv] += bc_mask.data()[i] * MLi * dm_du_n_fct / dt;
4168 for (
int i = 0; i < numDOFs; i++) {
4169 globalResidual.data()[offset_u + stride_u * i] += fluxCorrection.data()[i];
4171 for (
int i_n = 0; i_n < numDOFs_n; i_n++) {
4172 globalResidual.data()[offset_n + stride_n * i_n] += fluxCorrection_n.data()[i_n];
4191 std::vector<double> m_n_DOF(numDOFs_n, 0.0);
4192 std::vector<double> mn_n_DOF(numDOFs_n, 0.0);
4193 for (
int i_n = 0; i_n < numDOFs_n; i_n++) {
4194 const double sat = u_dof_n.data()[i_n];
4195 const double sat_old = u_dof_n_old.data()[i_n];
4196 m_n_DOF[i_n] = rho_n_phi_dof[i_n] * sat;
4197 mn_n_DOF[i_n] = rho_n_phi_dof_old[i_n] * sat_old;
4198 mn_n.data()[i_n] = mn_n_DOF[i_n];
4199 quantDOFs_n.data()[i_n] = 0.0;
4202 double diag_sumF = 0.0, diag_absF = 0.0;
4205 const bool have_gas_budget = (gas_budget_node.size() >= (
size_t)(6 * numDOFs_n));
4206 if (have_gas_budget)
4207 for (
int s = 0;
s < 6 * numDOFs_n; ++
s) gas_budget_node.data()[
s] = 0.0;
4209 std::vector<double> psi_n(numDOFs_n, 1.0);
4211 for (
int i_n = 0; i_n < numDOFs_n; i_n++) {
4212 const double zi = u_dof_n_old.data()[i_n];
4213 double num = 0.0, den = 0.0;
4214 for (
int offset = csrRowIndeces_n_DofLoops.data()[i_n];
4215 offset < csrRowIndeces_n_DofLoops.data()[i_n + 1]; offset++) {
4216 const int j_n = csrColumnOffsets_n_DofLoops.data()[offset];
4217 if (j_n == i_n)
continue;
4218 const double d = zi - u_dof_n_old.data()[j_n];
4222 const double alpha_i = fabs(
num) / (den + 1.0e-15);
4227 std::vector<int> node_iface(numDOFs_n, 0);
4229 std::vector<int> node_mat0(numDOFs_n, -1);
4230 for (
int eN = 0; eN < nElements_global; eN++) {
4231 const int mat_eN_if = elementMaterialTypes.data()[eN];
4232 for (
int a = 0; a < nDOF_trial_element; a++) {
4233 const int gN = u_l2g.data()[eN * nDOF_trial_element + a];
4234 if (node_mat0[gN] < 0) node_mat0[gN] = mat_eN_if;
4235 else if (node_mat0[gN] != mat_eN_if) node_iface[gN] = 1;
4239 auto comp1_offset = [&](
int i_n,
int j_n) ->
int {
4240 for (
int off = csrRowIndeces_n_DofLoops.data()[i_n];
4241 off < csrRowIndeces_n_DofLoops.data()[i_n + 1]; ++off)
4242 if (csrColumnOffsets_n_DofLoops.data()[off] == j_n)
return off;
4245 for (
int off = 0; off < NNZ_n; ++off) {
4246 dLow_n.data()[off] = 0.0;
4247 dEV_n.data()[off] = 0.0;
4248 dt_times_fH_minus_fL_n.data()[off] = 0.0;
4250 std::vector<double> node_Kdiag(numDOFs_n, 0.0);
4252 for (
int eN = 0; eN < nElements_global; eN++) {
4253 const int mat_eN = elementMaterialTypes.data()[eN];
4254 const double phi_eN = thetaR.data()[mat_eN] + thetaSR.data()[mat_eN];
4255 const double alpha_eN = alpha.data()[mat_eN];
4256 const double krn_end_eN = krn_end.data()[mat_eN];
4257 const double n_vg_eN =
n.data()[mat_eN];
4258 const double *KWs_eN = &KWs.data()[mat_eN *
nnz];
4262 double elementResidual_n[nDOF_test_element];
4267 double elementResidual_w[nDOF_test_element];
4268 double elementMass_n[nDOF_test_element];
4269 double u_n_local[nDOF_trial_element];
4270 double u_n_old_local[nDOF_trial_element];
4271 double elementJacobian_n_n[nDOF_test_element][nDOF_trial_element];
4272 double elementJacobian_n_w[nDOF_test_element][nDOF_trial_element];
4282 double elementJacobian_w_w[nDOF_test_element][nDOF_trial_element];
4283 double elementJacobian_w_n[nDOF_test_element][nDOF_trial_element];
4289 double elementTransport_n[nDOF_test_element][nDOF_trial_element];
4290 const int eN_nDOF_trial_element = eN * nDOF_trial_element;
4291 for (
int i = 0; i < nDOF_test_element; i++) {
4292 elementResidual_n[i] = 0.0;
4293 elementResidual_w[i] = 0.0;
4294 elementMass_n[i] = 0.0;
4295 for (
int j = 0; j < nDOF_trial_element; j++) {
4296 elementJacobian_n_n[i][j] = 0.0;
4297 elementJacobian_n_w[i][j] = 0.0;
4298 elementJacobian_w_w[i][j] = 0.0;
4299 elementJacobian_w_n[i][j] = 0.0;
4300 elementTransport_n[i][j] = 0.0;
4303 for (
int j = 0; j < nDOF_trial_element; j++) {
4304 u_n_local[j] = u_dof_n.data()[u_l2g_n.data()[eN_nDOF_trial_element + j]];
4305 u_n_old_local[j] = u_dof_n_old.data()[u_l2g_n.data()[eN_nDOF_trial_element + j]];
4307 for (
int k = 0; k < nQuadraturePoints_element; k++) {
4308 double jac[nSpace * nSpace], jacDet, jacInv[nSpace * nSpace], x_q, y_q, z_q;
4309 ck.calculateMapping_element(eN, k, mesh_dof.data(), mesh_l2g.data(),
4310 mesh_trial_ref.data(), mesh_grad_trial_ref.data(),
4311 jac, jacDet, jacInv, x_q, y_q, z_q);
4312 const double dV = std::fabs(jacDet) * dV_ref.data()[k];
4313 double u_grad_trial_qp[nDOF_trial_element * nSpace];
4314 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k * nDOF_trial_element * nSpace],
4315 jacInv, u_grad_trial_qp);
4326 for (
int i = 0; i < nDOF_test_element; i++) {
4327 const double test_i = u_test_ref.data()[k * nDOF_test_element + i];
4328 elementMass_n[i] += test_i * dV;
4329 for (
int j = 0; j < nDOF_trial_element; j++) {
4330 double K_trial_ij = 0.0;
4331 for (
int I = 0; I < nSpace; I++) {
4332 const double grad_Ni_I = u_grad_trial_qp[i * nSpace + I];
4333 for (
int ii = a_rowptr.data()[I]; ii < a_rowptr.data()[I + 1]; ii++) {
4334 const int J = a_colind.data()[ii];
4335 K_trial_ij += KWs_eN[ii] * u_grad_trial_qp[j * nSpace + J] * grad_Ni_I;
4338 elementTransport_n[i][j] += K_trial_ij * dV;
4350 for (
int i = 0; i < nDOF_test_element; i++) {
4351 const int gi = u_l2g_n.data()[eN * nDOF_test_element + i];
4352 node_Kdiag[gi] += fmax(0.0, elementTransport_n[i][i]);
4353 const double phiN_i = rho_n_phi_dof[gi];
4354 const double phiN_old_i= rho_n_phi_dof_old[gi];
4355 const double z_i = u_n_local[i];
4356 const double m_n_loc = phiN_i * z_i;
4357 const double m_n_old_loc = phiN_old_i * u_n_old_local[i];
4358 elementResidual_n[i] += elementMass_n[i] * (m_n_loc - m_n_old_loc) / dt;
4359 if (have_gas_budget)
4360 gas_budget_node.data()[0 * numDOFs_n + gi] += elementMass_n[i] * (m_n_loc - m_n_old_loc) / dt;
4362 elementJacobian_n_n[i][i] += elementMass_n[i] * (phiN_i + z_i * dphiN_dz_dof[gi]) / dt;
4364 elementJacobian_n_w[i][i] += elementMass_n[i] * (z_i * dphiN_dp_dof[gi]) / dt;
4365 if (inj_point_mode == 0) {
4366 const double Q_inj_n = injection_dof.data()[gi];
4367 elementResidual_n[i] -= elementMass_n[i] * Q_inj_n;
4368 if (have_gas_budget)
4369 gas_budget_node.data()[3 * numDOFs_n + gi] -= elementMass_n[i] * Q_inj_n;
4373 if (inj_point_mode == 1) {
4374 for (
int p = 0; p < inj_n_ports; p++) {
4375 if (inj_element.data()[p] == eN && inj_rate.data()[p] != 0.0) {
4376 const double qp = inj_rate.data()[p];
4377 for (
int i = 0; i < nDOF_test_element; i++) {
4378 const double w = inj_weight.data()[p * nDOF_test_element + i];
4379 elementResidual_n[i] -= qp *
w;
4380 if (have_gas_budget) {
4381 const int gii = u_l2g_n.data()[eN * nDOF_test_element + i];
4382 gas_budget_node.data()[3 * numDOFs_n + gii] -= qp *
w;
4390 const double S_wr_eN = thetaR.data()[mat_eN] / phi_eN;
4391 const double one_m_Sr_eN = 1.0 - S_wr_eN;
4392 const double Se_trap_L4415 = 1.0 - S_gr.data()[mat_eN] / one_m_Sr_eN;
4397 const double p_d_e = (alpha_eN > 0.0) ? (1.0 / alpha_eN) : 0.0;
4403 int gN_e[nDOF_trial_element];
4404 int zN_e[nDOF_trial_element];
4405 double uw_e[nDOF_trial_element], uw_old_e[nDOF_trial_element];
4406 double pc_e[nDOF_trial_element], pc_old_e[nDOF_trial_element];
4407 double dpc_dp_e[nDOF_trial_element], dpc_dz_e[nDOF_trial_element];
4408 double lam_g_e[nDOF_trial_element], lam_a_e[nDOF_trial_element];
4409 double lam_g_old_e[nDOF_trial_element], lam_a_old_e[nDOF_trial_element];
4410 double dlam_g_dp_e[nDOF_trial_element], dlam_g_dz_e[nDOF_trial_element];
4411 double dlam_a_dp_e[nDOF_trial_element], dlam_a_dz_e[nDOF_trial_element];
4416 double lwg_e[nDOF_trial_element], lwa_e[nDOF_trial_element];
4417 double lwg_old_e[nDOF_trial_element], lwa_old_e[nDOF_trial_element];
4418 double dlwg_dp_e[nDOF_trial_element], dlwg_dz_e[nDOF_trial_element];
4419 double dlwa_dp_e[nDOF_trial_element], dlwa_dz_e[nDOF_trial_element];
4420 double rgm_e[nDOF_trial_element], ram_e[nDOF_trial_element];
4421 double rgm_old_e[nDOF_trial_element], ram_old_e[nDOF_trial_element];
4422 double drgm_dp_e[nDOF_trial_element], drgm_dz_e[nDOF_trial_element];
4423 double dram_dp_e[nDOF_trial_element], dram_dz_e[nDOF_trial_element];
4424 double Sg_e[nDOF_trial_element], Sg_old_e[nDOF_trial_element];
4425 double dSg_dp_e[nDOF_trial_element], dSg_dz_e[nDOF_trial_element];
4426 double dlg_dz_o[nDOF_trial_element], dla_dz_o[nDOF_trial_element], dpc_dz_o[nDOF_trial_element];
4427 const double cg_eN = krn_end_eN / mu_n;
4429 for (
int a = 0; a < nDOF_trial_element; a++) {
4430 const int gN = u_l2g.data()[eN_nDOF_trial_element + a];
4431 const int zN = u_l2g_n.data()[eN_nDOF_trial_element + a];
4434 const double p_a = u_dof.data()[gN];
4435 const double z_a = u_dof_n.data()[zN];
4436 const double p_a_o = u_dof_old.data()[gN];
4437 const double z_a_o = u_dof_n_old.data()[zN];
4438 uw_e[a] = p_a; uw_old_e[a] = p_a_o;
4440 const double z_cl = fmin(fmax(z_a, 1.0e-8), 1.0 - 1.0e-8);
4441 const double p_cl = fmax(p_a, 1.0e2);
4444 Sg_e[a] =
f.S_g; dSg_dp_e[a] =
f.dS_g_dp; dSg_dz_e[a] =
f.dS_g_dz;
4445 const double Sa = 1.0 -
f.S_g;
4446 const double Se_raw = (Sa - S_wr_eN) / one_m_Sr_eN;
4447 double Se, dSe_dp, dSe_dz;
4448 if (Se_raw <= 0.0) { Se = 0.0; dSe_dp = 0.0; dSe_dz = 0.0; }
4449 else if (Se_raw >= 1.0) { Se = 1.0; dSe_dp = 0.0; dSe_dz = 0.0; }
4450 else { Se = Se_raw; dSe_dp = -
f.dS_g_dp/one_m_Sr_eN; dSe_dz = -
f.dS_g_dz/one_m_Sr_eN; }
4451 double krn=0,dkrn=0,krw=0,dkrw=0,thW=0,DthW=0,pc=0,dpc_dSe=0,d2pc=0;
4461 pc_e[a] = pc; dpc_dp_e[a] = dpc_dSe*dSe_dp; dpc_dz_e[a] = dpc_dSe*dSe_dz;
4464 rgm_e[a] =
f.rho_g*Mg; ram_e[a] =
f.rho_a*Ma;
4465 drgm_dp_e[a] =
f.drho_g_dp*Mg +
f.rho_g*
f.dY_dp*dMm_e; drgm_dz_e[a] =
f.rho_g*
f.dY_dz*dMm_e;
4466 dram_dp_e[a] =
f.drho_a_dp*Ma +
f.rho_a*
f.dX_dp*dMm_e; dram_dz_e[a] =
f.drho_a_dz*Ma +
f.rho_a*
f.dX_dz*dMm_e;
4468 lam_g_e[a] = cg_eN*
f.rho_g*
f.Y*krn;
4469 dlam_g_dp_e[a] = cg_eN*(
f.drho_g_dp*
f.Y*krn +
f.rho_g*
f.dY_dp*krn +
f.rho_g*
f.Y*dkrn*dSe_dp);
4470 dlam_g_dz_e[a] = cg_eN*(
f.rho_g*
f.dY_dz*krn +
f.rho_g*
f.Y*dkrn*dSe_dz);
4471 lam_a_e[a] =
f.rho_a*
f.X*krw;
4472 dlam_a_dp_e[a] =
f.drho_a_dp*
f.X*krw +
f.rho_a*
f.dX_dp*krw +
f.rho_a*
f.X*dkrw*dSe_dp;
4473 dlam_a_dz_e[a] =
f.drho_a_dz*
f.X*krw +
f.rho_a*
f.dX_dz*krw +
f.rho_a*
f.X*dkrw*dSe_dz;
4476 const double yw=1.0-
f.Y, dyw_dp=-
f.dY_dp, dyw_dz=-
f.dY_dz;
4477 const double xw=1.0-
f.X, dxw_dp=-
f.dX_dp, dxw_dz=-
f.dX_dz;
4478 lwg_e[a] = cg_eN*
f.rho_g*yw*krn;
4479 dlwg_dp_e[a] = cg_eN*(
f.drho_g_dp*yw*krn +
f.rho_g*dyw_dp*krn +
f.rho_g*yw*dkrn*dSe_dp);
4480 dlwg_dz_e[a] = cg_eN*(
f.rho_g*dyw_dz*krn +
f.rho_g*yw*dkrn*dSe_dz);
4481 lwa_e[a] =
f.rho_a*xw*krw;
4482 dlwa_dp_e[a] =
f.drho_a_dp*xw*krw +
f.rho_a*dxw_dp*krw +
f.rho_a*xw*dkrw*dSe_dp;
4483 dlwa_dz_e[a] =
f.drho_a_dz*xw*krw +
f.rho_a*dxw_dz*krw +
f.rho_a*xw*dkrw*dSe_dz;
4486 const double z_cl_o = fmin(fmax(z_a_o, 1.0e-8), 1.0 - 1.0e-8);
4487 const double p_cl_o = fmax(p_a_o, 1.0e2);
4490 Sg_old_e[a] = fo.
S_g;
4491 const double Sa_o = 1.0 - fo.
S_g;
4492 const double Se_o_raw = (Sa_o - S_wr_eN) / one_m_Sr_eN;
4493 const double Se_o = Se_o_raw <= 0.0 ? 0.0 : (Se_o_raw >= 1.0 ? 1.0 : Se_o_raw);
4494 double krn_o=0,dkrn_o=0,krw_o=0,dkrw_o=0,thW_o=0,DthW_o=0,pc_o=0,dpc_o=0,d2pc_o=0;
4507 rgm_old_e[a] = fo.
rho_g*Mg_o; ram_old_e[a] = fo.
rho_a*Ma_o;
4508 lam_g_old_e[a] = cg_eN*fo.
rho_g*fo.
Y*krn_o;
4509 lam_a_old_e[a] = fo.
rho_a*fo.
X*krw_o;
4510 lwg_old_e[a] = cg_eN*fo.
rho_g*(1.0-fo.
Y)*krn_o;
4511 lwa_old_e[a] = fo.
rho_a*(1.0-fo.
X)*krw_o;
4512 const double dSe_o = (Se_o_raw>0.0 && Se_o_raw<1.0) ? -fo.
dS_g_dz/one_m_Sr_eN : 0.0;
4513 dlg_dz_o[a] = cg_eN*(fo.
rho_g*fo.
dY_dz*krn_o + fo.
rho_g*fo.
Y*dkrn_o*dSe_o);
4515 dpc_dz_o[a] = dpc_o*dSe_o;
4521 for (
int i = 0; i < nDOF_test_element; i++) {
4522 for (
int j = 0; j < nDOF_trial_element; j++) {
4523 if (i == j)
continue;
4524 const double tau = fmax(0.0, -elementTransport_n[i][j]);
4525 if (tau == 0.0)
continue;
4526 double g_dot_dx = 0.0;
4527 for (
int I = 0; I < nSpace; I++) {
4528 g_dot_dx += gravity.data()[I] * (mesh_dof.data()[gN_e[j] * 3 + I]
4529 - mesh_dof.data()[gN_e[i] * 3 + I]);
4533 const double rgm_edge = 0.5 * (rgm_e[i] + rgm_e[j]);
4534 const double rgm_edge_old = 0.5 * (rgm_old_e[i] + rgm_old_e[j]);
4535 const double dPhi_g = (uw_e[j] + pc_e[j]) - (uw_e[i] + pc_e[i])
4536 - rgm_edge * g_dot_dx;
4537 const double dPhi_g_old = (uw_old_e[j] + pc_old_e[j]) - (uw_old_e[i] + pc_old_e[i])
4538 - rgm_edge_old * g_dot_dx;
4539 const bool up_i_g = (dPhi_g <= 0.0);
4540 const bool up_i_g_old = (dPhi_g_old <= 0.0);
4541 const double lam_g_up = up_i_g ? lam_g_e[i] : lam_g_e[j];
4542 const double lam_g_up_old = up_i_g_old ? lam_g_old_e[i] : lam_g_old_e[j];
4550 const double Fg = Theta * tau * lam_g_up * dPhi_g
4551 + (1.0 - Theta) * tau * lam_g_up_old * dPhi_g_old;
4554 const double ram_edge = 0.5 * (ram_e[i] + ram_e[j]);
4555 const double ram_edge_old = 0.5 * (ram_old_e[i] + ram_old_e[j]);
4556 const double dPhi_a = uw_e[j] - uw_e[i] - ram_edge * g_dot_dx;
4557 const double dPhi_a_old = uw_old_e[j] - uw_old_e[i] - ram_edge_old * g_dot_dx;
4558 const bool up_i_a = (dPhi_a <= 0.0);
4559 const bool up_i_a_old = (dPhi_a_old <= 0.0);
4560 const double lam_a_up = up_i_a ? lam_a_e[i] : lam_a_e[j];
4561 const double lam_a_up_old = up_i_a_old ? lam_a_old_e[i] : lam_a_old_e[j];
4562 const double Fa = Theta * tau * lam_a_up * dPhi_a
4563 + (1.0 - Theta) * tau * lam_a_up_old * dPhi_a_old;
4565 elementResidual_n[i] -= Fg + Fa;
4566 if (have_gas_budget)
4567 gas_budget_node.data()[1 * numDOFs_n + zN_e[i]] -= Fg + Fa;
4568 diag_sumF += Fg + Fa;
4569 diag_absF += std::fabs(Fg + Fa);
4571 const double ddPhig_dpi = -1.0 - dpc_dp_e[i] - 0.5 * drgm_dp_e[i] * g_dot_dx;
4572 const double ddPhig_dzi = - dpc_dz_e[i] - 0.5 * drgm_dz_e[i] * g_dot_dx;
4573 const double ddPhig_dpj = +1.0 + dpc_dp_e[j] - 0.5 * drgm_dp_e[j] * g_dot_dx;
4574 const double ddPhig_dzj = + dpc_dz_e[j] - 0.5 * drgm_dz_e[j] * g_dot_dx;
4575 const double ddPhia_dpi = -1.0 - 0.5 * dram_dp_e[i] * g_dot_dx;
4576 const double ddPhia_dzi = - 0.5 * dram_dz_e[i] * g_dot_dx;
4577 const double ddPhia_dpj = +1.0 - 0.5 * dram_dp_e[j] * g_dot_dx;
4578 const double ddPhia_dzj = - 0.5 * dram_dz_e[j] * g_dot_dx;
4579 const double Tt = Theta * tau;
4580 double dF_dpi = 0.0, dF_dzi = 0.0, dF_dpj = 0.0, dF_dzj = 0.0;
4582 dF_dpi += Tt*lam_g_up*ddPhig_dpi; dF_dzi += Tt*lam_g_up*ddPhig_dzi;
4583 dF_dpj += Tt*lam_g_up*ddPhig_dpj; dF_dzj += Tt*lam_g_up*ddPhig_dzj;
4585 if (up_i_g) { dF_dpi += Tt*dlam_g_dp_e[i]*dPhi_g; dF_dzi += Tt*dlam_g_dz_e[i]*dPhi_g; }
4586 else { dF_dpj += Tt*dlam_g_dp_e[j]*dPhi_g; dF_dzj += Tt*dlam_g_dz_e[j]*dPhi_g; }
4589 dF_dpi += Tt*lam_a_up*ddPhia_dpi; dF_dzi += Tt*lam_a_up*ddPhia_dzi;
4590 dF_dpj += Tt*lam_a_up*ddPhia_dpj; dF_dzj += Tt*lam_a_up*ddPhia_dzj;
4592 if (up_i_a) { dF_dpi += Tt*dlam_a_dp_e[i]*dPhi_a; dF_dzi += Tt*dlam_a_dz_e[i]*dPhi_a; }
4593 else { dF_dpj += Tt*dlam_a_dp_e[j]*dPhi_a; dF_dzj += Tt*dlam_a_dz_e[j]*dPhi_a; }
4595 elementJacobian_n_n[i][i] += -dF_dzi;
4596 elementJacobian_n_n[i][j] += -dF_dzj;
4597 elementJacobian_n_w[i][i] += -dF_dpi;
4598 elementJacobian_n_w[i][j] += -dF_dpj;
4609 const double lwg_up = up_i_g ? lwg_e[i] : lwg_e[j];
4610 const double lwg_up_old = up_i_g_old ? lwg_old_e[i] : lwg_old_e[j];
4611 const double lwa_up = up_i_a ? lwa_e[i] : lwa_e[j];
4612 const double lwa_up_old = up_i_a_old ? lwa_old_e[i] : lwa_old_e[j];
4613 const double Fwg = Theta * tau * lwg_up * dPhi_g
4614 + (1.0 - Theta) * tau * lwg_up_old * dPhi_g_old;
4615 const double Fwa = Theta * tau * lwa_up * dPhi_a
4616 + (1.0 - Theta) * tau * lwa_up_old * dPhi_a_old;
4617 elementResidual_w[i] -= Fwg + Fwa;
4619 double dFw_dpi=0.0, dFw_dzi=0.0, dFw_dpj=0.0, dFw_dzj=0.0;
4620 dFw_dpi += Tt*lwg_up*ddPhig_dpi; dFw_dzi += Tt*lwg_up*ddPhig_dzi;
4621 dFw_dpj += Tt*lwg_up*ddPhig_dpj; dFw_dzj += Tt*lwg_up*ddPhig_dzj;
4622 if (up_i_g) { dFw_dpi += Tt*dlwg_dp_e[i]*dPhi_g; dFw_dzi += Tt*dlwg_dz_e[i]*dPhi_g; }
4623 else { dFw_dpj += Tt*dlwg_dp_e[j]*dPhi_g; dFw_dzj += Tt*dlwg_dz_e[j]*dPhi_g; }
4624 dFw_dpi += Tt*lwa_up*ddPhia_dpi; dFw_dzi += Tt*lwa_up*ddPhia_dzi;
4625 dFw_dpj += Tt*lwa_up*ddPhia_dpj; dFw_dzj += Tt*lwa_up*ddPhia_dzj;
4626 if (up_i_a) { dFw_dpi += Tt*dlwa_dp_e[i]*dPhi_a; dFw_dzi += Tt*dlwa_dz_e[i]*dPhi_a; }
4627 else { dFw_dpj += Tt*dlwa_dp_e[j]*dPhi_a; dFw_dzj += Tt*dlwa_dz_e[j]*dPhi_a; }
4634 const int fi = r_l2g.data()[eN_nDOF_trial_element + i];
4635 const double mwi = bc_mask.data()[fi];
4636 elementJacobian_w_w[i][i] += mwi * (-dFw_dpi);
4637 elementJacobian_w_w[i][j] += mwi * (-dFw_dpj);
4638 elementJacobian_w_n[i][i] += mwi * (-dFw_dzi);
4639 elementJacobian_w_n[i][j] += mwi * (-dFw_dzj);
4657 const double si_g = std::fabs(dlg_dz_o[i])*std::fabs(dPhi_g_old)
4658 + std::fabs(lam_g_old_e[i])*std::fabs(dpc_dz_o[i]);
4659 const double si_a = std::fabs(dla_dz_o[i])*std::fabs(dPhi_a_old);
4660 const double sj_g = std::fabs(dlg_dz_o[j])*std::fabs(dPhi_g_old)
4661 + std::fabs(lam_g_old_e[j])*std::fabs(dpc_dz_o[j]);
4662 const double sj_a = std::fabs(dla_dz_o[j])*std::fabs(dPhi_a_old);
4663 const double psi_edge = fmax(psi_n[zN_e[i]], psi_n[zN_e[j]]);
4664 const double dEV = tau * (
cE*psi_edge*fmax(si_g,sj_g)
4665 +
cE*psi_edge*fmax(si_a,sj_a) );
4666 const double dLow = tau * ( fmax(si_g,sj_g) + fmax(si_a,sj_a) );
4667 const double dHi = fmin(dEV, dLow);
4668 const double dResid = (FCT_n == 1) ? dLow : dEV;
4670 const double Fv = dResid*(u_dof_n.data()[zN_e[j]] - u_dof_n.data()[zN_e[i]]);
4671 elementResidual_n[i] -= Fv;
4672 elementJacobian_n_n[i][i] += dResid;
4673 elementJacobian_n_n[i][j] -= dResid;
4676 const int off_n = comp1_offset(zN_e[i], zN_e[j]);
4678 const double dz = u_dof_n.data()[zN_e[j]] - u_dof_n.data()[zN_e[i]];
4679 dLow_n.data()[off_n] += dLow;
4680 dEV_n.data()[off_n] += dHi;
4681 dt_times_fH_minus_fL_n.data()[off_n] += dt * (dLow - dHi) * dz;
4689 for (
int i = 0; i < nDOF_test_element; i++) {
4690 const int eN_i = eN * nDOF_test_element + i;
4691 const int gi = u_l2g_n.data()[eN_i];
4692 globalResidual.data()[offset_n + stride_n * gi] += elementResidual_n[i];
4693 if (have_gas_budget)
4694 gas_budget_node.data()[5 * numDOFs_n + gi] += elementResidual_n[i];
4696 const int fi_w = r_l2g.data()[eN_nDOF_trial_element + i];
4697 globalResidual.data()[offset_u + stride_u * fi_w] += bc_mask.data()[fi_w] * elementResidual_w[i];
4698 for (
int j = 0; j < nDOF_trial_element; j++) {
4699 const int eN_i_j = eN_i * nDOF_trial_element + j;
4700 globalJacobian.data()[csrRowIndeces_n_n.data()[eN_i] + csrColumnOffsets_n_n.data()[eN_i_j]]
4701 += elementJacobian_n_n[i][j];
4702 globalJacobian.data()[csrRowIndeces_n_w.data()[eN_i] + csrColumnOffsets_n_w.data()[eN_i_j]]
4703 += elementJacobian_n_w[i][j];
4704 globalJacobian.data()[csrRowIndeces_w_w.data()[eN_i] + csrColumnOffsets_w_w.data()[eN_i_j]]
4705 += elementJacobian_w_w[i][j];
4706 globalJacobian.data()[csrRowIndeces_w_n.data()[eN_i] + csrColumnOffsets_w_n.data()[eN_i_j]]
4707 += elementJacobian_w_n[i][j];
4712 for (
int i_n = 0; i_n < numDOFs_n; i_n++) {
4713 mLow_n.data()[i_n] = m_n_DOF[i_n];
4714 mDotLow_n.data()[i_n] = (m_n_DOF[i_n] - mn_n.data()[i_n]) / dt;
4718 const double mu_n_loc = mu_n;
4719 int n_drop_offdiag = 0, n_drop_pcol = 0;
4720 for (
int ip = 0; ip < n_interface_pairs; ip++) {
4721 const int nodeN = interface_pairs.data()[5 * ip + 0];
4722 const int zc[2] = { interface_pairs.data()[5 * ip + 1], interface_pairs.data()[5 * ip + 3] };
4723 const int mc[2] = { interface_pairs.data()[5 * ip + 2], interface_pairs.data()[5 * ip + 4] };
4724 const double p_node = u_dof.data()[nodeN];
4725 double pc2[2], dpc_dz2[2], lam_g2[2], dlam_g_dz2[2], C2[2], dC_dz2[2];
4726 double dpc_dp2[2], dlam_g_dp2[2], dC_dp2[2];
4727 for (
int s = 0;
s < 2;
s++) {
4728 const int mat = mc[
s];
4729 const double phi_m = thetaR.data()[mat] + thetaSR.data()[mat];
4730 const double S_wr_m = thetaR.data()[mat] / phi_m;
4731 const double oneSr_m = 1.0 - S_wr_m;
4732 const double alpha_m = alpha.data()[mat];
4733 const double n_m =
n.data()[mat];
4734 const double cg_m = krn_end.data()[mat] / mu_n_loc;
4735 const double Setrap_m= 1.0 - S_gr.data()[mat] / oneSr_m;
4736 const double z_cl = fmin(fmax(u_dof_n.data()[zc[
s]], 1.0e-8), 1.0 - 1.0e-8);
4737 const double p_cl = fmax(p_node, 1.0e2);
4740 const double Sa = 1.0 -
f.S_g;
4741 const double Se_raw = (Sa - S_wr_m) / oneSr_m;
4742 double Se, dSe_dz, dSe_dp;
4743 if (Se_raw <= 0.0) { Se = 0.0; dSe_dz = 0.0; dSe_dp = 0.0; }
4744 else if (Se_raw >= 1.0) { Se = 1.0; dSe_dz = 0.0; dSe_dp = 0.0; }
4745 else { Se = Se_raw; dSe_dz = -
f.dS_g_dz / oneSr_m; dSe_dp = -
f.dS_g_dp / oneSr_m; }
4746 double krn = 0, dkrn = 0, kpc = 0, dpc_dSe = 0, d2pc = 0, thw = 0, Dthw = 0, krw = 0, dkrw = 0;
4755 dpc_dz2[
s] = dpc_dSe * dSe_dz;
4756 dpc_dp2[
s] = dpc_dSe * dSe_dp;
4757 lam_g2[
s] = cg_m *
f.rho_g *
f.Y * krn;
4758 dlam_g_dz2[
s]= cg_m * (
f.rho_g *
f.dY_dz * krn +
f.rho_g *
f.Y * dkrn * dSe_dz);
4761 dlam_g_dp2[
s]= cg_m * (
f.drho_g_dp *
f.Y * krn +
f.rho_g *
f.dY_dp * krn
4762 +
f.rho_g *
f.Y * dkrn * dSe_dp);
4763 C2[
s] =
f.rho_a *
f.X;
4764 dC_dz2[
s] =
f.drho_a_dz *
f.X +
f.rho_a *
f.dX_dz;
4765 dC_dp2[
s] =
f.drho_a_dp *
f.X +
f.rho_a *
f.dX_dp;
4767 const int z_a = zc[0], z_b = zc[1];
4768 const double Ka = node_Kdiag[z_a], Kb = node_Kdiag[z_b];
4769 const double tau_if = (Ka + Kb > 0.0) ? (2.0 * Ka * Kb / (Ka + Kb)) : 0.0;
4770 if (tau_if == 0.0)
continue;
4771 const double dPhi_g = pc2[0] - pc2[1];
4772 const bool up_a = (dPhi_g >= 0.0);
4773 const double lam_up = up_a ? lam_g2[0] : lam_g2[1];
4774 const double dlam_up_dza = up_a ? dlam_g_dz2[0] : 0.0;
4775 const double dlam_up_dzb = up_a ? 0.0 : dlam_g_dz2[1];
4776 const double Fg = tau_if * lam_up * dPhi_g;
4777 const double Fd = D_m * tau_if * (C2[0] - C2[1]);
4778 const double F = Fg + Fd;
4779 const double dF_dza = tau_if * (dlam_up_dza * dPhi_g + lam_up * dpc_dz2[0]) + D_m * tau_if * ( dC_dz2[0]);
4780 const double dF_dzb = tau_if * (dlam_up_dzb * dPhi_g - lam_up * dpc_dz2[1]) + D_m * tau_if * (-dC_dz2[1]);
4785 const double dlam_up_dp = up_a ? dlam_g_dp2[0] : dlam_g_dp2[1];
4786 const double dF_dp = tau_if * (dlam_up_dp * dPhi_g + lam_up * (dpc_dp2[0] - dpc_dp2[1]))
4787 + D_m * tau_if * (dC_dp2[0] - dC_dp2[1]);
4789 globalResidual.data()[offset_n + stride_n * z_a] += F;
4790 globalResidual.data()[offset_n + stride_n * z_b] -= F;
4793 const int caa = comp1_offset(z_a, z_a), cbb = comp1_offset(z_b, z_b);
4794 if (caa >= 0) globalJacobian.data()[comp1_full_offsets.data()[caa]] += dF_dza;
4795 if (cbb >= 0) globalJacobian.data()[comp1_full_offsets.data()[cbb]] += -dF_dzb;
4799 const int cab = comp1_iface_offsets.data()[2 * ip + 0];
4800 const int cba = comp1_iface_offsets.data()[2 * ip + 1];
4801 if (cab >= 0) globalJacobian.data()[cab] += dF_dzb;
4802 else ++n_drop_offdiag;
4803 if (cba >= 0) globalJacobian.data()[cba] += -dF_dza;
4804 else ++n_drop_offdiag;
4807 const int oa = comp10_full_offsets.data()[2 * ip + 0];
4808 const int ob = comp10_full_offsets.data()[2 * ip + 1];
4809 if (oa >= 0) globalJacobian.data()[oa] += dF_dp;
4811 if (ob >= 0) globalJacobian.data()[ob] += -dF_dp;
4814 if (n_drop_offdiag > 0 || n_drop_pcol > 0) {
4815 std::cerr <<
"[m_comp_co2 split_z] WARNING: dropped interface tangents -- "
4816 << n_drop_offdiag <<
" z_a<->z_b off-diagonal, "
4817 << n_drop_pcol <<
" z<->p_node (1,0) slot(s) absent "
4818 <<
"(stale/incomplete Jacobian sparsity -> degraded Newton)"
4822 if (split_anchor_alpha > 0.0 && dt > 0.0) {
4823 const double Sg_tol = split_anchor_Sg_tol;
4824 const double X_tol = split_anchor_X_tol;
4828 if (split_anchor_layer1 != 0)
4829 for (
int i = 0; i < numDOFs_n; i++) {
4830 if (ML_n[i] <= 0.0)
continue;
4831 if (Sg_dof_old[i] >= Sg_tol || X_dof_old[i] >= X_tol)
continue;
4832 const double cap_i = rho_n_phi_dof_old[i] * ML_n[i];
4833 const double zi = u_dof_n.data()[i];
4835 for (
int off = csrRowIndeces_n_DofLoops.data()[i];
4836 off < csrRowIndeces_n_DofLoops.data()[i + 1]; ++off) {
4837 const int j = csrColumnOffsets_n_DofLoops.data()[off];
4838 if (j == i || ML_n[j] <= 0.0)
continue;
4839 if (Sg_dof_old[j] >= Sg_tol || X_dof_old[j] >= X_tol)
continue;
4840 const double cap_j = rho_n_phi_dof_old[j] * ML_n[j];
4841 const double lam_ij = split_anchor_alpha * fmin(cap_i, cap_j) / dt;
4842 if (lam_ij <= 0.0)
continue;
4844 globalResidual.data()[offset_n + stride_n * i] += lam_ij * (zi - u_dof_n.data()[j]);
4847 globalJacobian.data()[comp1_full_offsets.data()[off]] += -lam_ij;
4851 const int ii = comp1_offset(i, i);
4852 if (ii >= 0) globalJacobian.data()[comp1_full_offsets.data()[ii]] += diag;
4856 for (
int ip = 0; ip < n_interface_pairs; ip++) {
4857 const int z_a = interface_pairs.data()[5 * ip + 1];
4858 const int z_b = interface_pairs.data()[5 * ip + 3];
4859 if (ML_n[z_a] <= 0.0 || ML_n[z_b] <= 0.0)
continue;
4862 if (Sg_dof_old[z_a] >= Sg_tol || X_dof_old[z_a] >= X_tol)
continue;
4863 if (Sg_dof_old[z_b] >= Sg_tol || X_dof_old[z_b] >= X_tol)
continue;
4864 const double cap_a = rho_n_phi_dof_old[z_a] * ML_n[z_a];
4865 const double cap_b = rho_n_phi_dof_old[z_b] * ML_n[z_b];
4866 const double lam_s = split_anchor_alpha * fmin(cap_a, cap_b) / dt;
4867 if (lam_s <= 0.0)
continue;
4868 const double F_s = lam_s * (u_dof_n.data()[z_a] - u_dof_n.data()[z_b]);
4869 globalResidual.data()[offset_n + stride_n * z_a] += F_s;
4870 globalResidual.data()[offset_n + stride_n * z_b] -= F_s;
4872 const int caa = comp1_offset(z_a, z_a), cbb = comp1_offset(z_b, z_b);
4873 if (caa >= 0) globalJacobian.data()[comp1_full_offsets.data()[caa]] += lam_s;
4874 if (cbb >= 0) globalJacobian.data()[comp1_full_offsets.data()[cbb]] += lam_s;
4877 const int cab = comp1_iface_offsets.data()[2 * ip + 0];
4878 const int cba = comp1_iface_offsets.data()[2 * ip + 1];
4879 if (cab >= 0) globalJacobian.data()[cab] += -lam_s;
4880 if (cba >= 0) globalJacobian.data()[cba] += -lam_s;
4885 if (gas_diag.size() >= 4) {
4886 gas_diag.data()[0] = 0.0;
4887 gas_diag.data()[1] = 0.0;
4888 gas_diag.data()[2] = diag_sumF;
4889 gas_diag.data()[3] = diag_absF;
4926 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++) {
4927 const int ebN = exteriorElementBoundariesArray.data()[ebNE];
4928 const int eN = elementBoundaryElementsArray.data()[ebN * 2 + 0];
4929 const int ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN * 2 + 0];
4930 const int eN_nDOF_trial_element = eN * nDOF_trial_element;
4931 const int mat_eN = elementMaterialTypes.data()[eN];
4932 const double phi_eN = thetaR.data()[mat_eN] + thetaSR.data()[mat_eN];
4933 const double alpha_eN = alpha.data()[mat_eN];
4934 const double krn_end_eN = krn_end.data()[mat_eN];
4935 const double n_vg_eN =
n.data()[mat_eN];
4936 const double *KWs_eN = &KWs.data()[mat_eN *
nnz];
4937 const double S_wr_loc = thetaR.data()[mat_eN] / phi_eN;
4938 const double one_m_Sr_loc = 1.0 - S_wr_loc;
4939 const double Se_trap_L4846 = 1.0 - S_gr.data()[mat_eN] / one_m_Sr_loc;
4941 double elementResidual_n_eb[nDOF_test_element];
4942 double elementJacobian_n_n_eb[nDOF_test_element][nDOF_trial_element];
4943 double elementJacobian_n_w_eb[nDOF_test_element][nDOF_trial_element];
4944 for (
int i = 0; i < nDOF_test_element; i++) {
4945 elementResidual_n_eb[i] = 0.0;
4946 for (
int j = 0; j < nDOF_trial_element; j++) {
4947 elementJacobian_n_n_eb[i][j] = 0.0;
4948 elementJacobian_n_w_eb[i][j] = 0.0;
4952 for (
int kb = 0; kb < nQuadraturePoints_elementBoundary; kb++) {
4953 const int ebNE_kb = ebNE * nQuadraturePoints_elementBoundary + kb;
4954 const int ebN_local_kb = ebN_local * nQuadraturePoints_elementBoundary + kb;
4955 const int ebN_local_kb_nSpace = ebN_local_kb * nSpace;
4957 double jac_ext[nSpace * nSpace], jacDet_ext, jacInv_ext[nSpace * nSpace];
4958 double boundaryJac_b[nSpace * (nSpace - 1)];
4959 double metricTensor_b[(nSpace - 1) * (nSpace - 1)];
4960 double metricTensorDetSqrt_b, dS_eb, normal_b[3];
4961 double xt_b, yt_b, zt_b, integralScaling_b;
4962 double x_eb, y_eb, z_eb;
4963 ck.calculateMapping_elementBoundary(eN, ebN_local, kb, ebN_local_kb,
4964 mesh_dof.data(), mesh_l2g.data(), mesh_trial_trace_ref.data(),
4965 mesh_grad_trial_trace_ref.data(), boundaryJac_ref.data(),
4966 jac_ext, jacDet_ext, jacInv_ext, boundaryJac_b, metricTensor_b,
4967 metricTensorDetSqrt_b, normal_ref.data(), normal_b,
4969 ck.calculateMappingVelocity_elementBoundary(eN, ebN_local, kb, ebN_local_kb,
4970 mesh_velocity_dof.data(), mesh_l2g.data(), mesh_trial_trace_ref.data(),
4971 xt_b, yt_b, zt_b, normal_b, boundaryJac_b, metricTensor_b,
4973 dS_eb = ((1.0 - MOVING_DOMAIN) * metricTensorDetSqrt_b
4974 + MOVING_DOMAIN * integralScaling_b) * dS_ref.data()[kb];
4977 double u_grad_trial_trace_b[nDOF_trial_element * nSpace];
4978 ck.gradTrialFromRef(
4979 &u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace * nDOF_trial_element],
4980 jacInv_ext, u_grad_trial_trace_b);
4981 double u_w_ext_b = 0.0, u_n_ext_b = 0.0;
4982 double grad_u_w_ext_b[nSpace], grad_u_n_ext_b[nSpace];
4983 ck.valFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element],
4984 &u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element],
4986 ck.valFromDOF(u_dof_n.data(), &u_l2g.data()[eN_nDOF_trial_element],
4987 &u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element],
4989 ck.gradFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element],
4990 u_grad_trial_trace_b, grad_u_w_ext_b);
4991 ck.gradFromDOF(u_dof_n.data(), &u_l2g.data()[eN_nDOF_trial_element],
4992 u_grad_trial_trace_b, grad_u_n_ext_b);
4995 const int isDir_n = isDOFBoundary_n.data()[ebNE_kb];
4996 const double bc_u_n_ext_b = isDir_n * ebqe_bc_u_n_ext.data()[ebNE_kb]
4997 + (1 - isDir_n) * u_n_ext_b;
5000 const double z_clb = fmin(fmax(u_n_ext_b, 1.0e-8), 1.0 - 1.0e-8);
5001 const double p_clb = fmax(u_w_ext_b, 1.0e2);
5004 const double S_g_b = fsb.
S_g, Sa_b = 1.0 - S_g_b;
5006 const double Se_raw_b = (Sa_b - S_wr_loc) / one_m_Sr_loc;
5007 double Se_b, dSe_dp_b, dSe_dz_b;
5008 if (Se_raw_b <= 0.0) { Se_b = 0.0; dSe_dp_b = 0.0; dSe_dz_b = 0.0; }
5009 else if (Se_raw_b >= 1.0) { Se_b = 1.0; dSe_dp_b = 0.0; dSe_dz_b = 0.0; }
5010 else { Se_b = Se_raw_b; dSe_dp_b = -fsb.
dS_g_dp/one_m_Sr_loc; dSe_dz_b = -fsb.
dS_g_dz/one_m_Sr_loc; }
5011 double KWr_b=0.0, DKWr_b=0.0, thW_b=0.0, DthW_b=0.0, KNr_b=0.0, DKNr_b=0.0;
5012 double pc_b=0.0, dpc_dSe_b=0.0, d2pc_b=0.0;
5015 thetaR.data()[mat_eN], thetaSR.data()[mat_eN], thW_b, DthW_b, KWr_b, DKWr_b);
5020 thetaR.data()[mat_eN], thetaSR.data()[mat_eN], thW_b, DthW_b, KWr_b, DKWr_b);
5024 KNr_b *= krn_end_eN; DKNr_b *= krn_end_eN;
5025 const double pcp_b = dpc_dSe_b / one_m_Sr_loc;
5026 const double dpcp_dp_b = (d2pc_b / one_m_Sr_loc) * dSe_dp_b;
5027 const double dpcp_dz_b = (d2pc_b / one_m_Sr_loc) * dSe_dz_b;
5032 const double rho_g_mass_b = fsb.
rho_g*Mbar_g_b, rho_a_mass_b = fsb.
rho_a*Mbar_a_b;
5034 const double drgm_dz_b = fsb.
rho_g*fsb.
dY_dz*dMm_b;
5038 const double Ag = fsb.
rho_g*fsb.
Y, Aa = fsb.
rho_a*fsb.
X;
5045 double ug_b[nSpace], ua_b[nSpace];
5046 double dug_dp_b[nSpace], dug_dz_b[nSpace], dua_dp_b[nSpace], dua_dz_b[nSpace];
5047 for (
int I = 0; I < nSpace; I++) {
5048 double ugI=0.0, uaI=0.0, dugp=0.0, dugz=0.0, duap=0.0, duaz=0.0;
5049 for (
int ii = a_rowptr.data()[I]; ii < a_rowptr.data()[I + 1]; ii++) {
5050 const int J = a_colind.data()[ii];
5051 const double Kii = KWs_eN[ii];
5052 const double Mob_g = KNr_b*Kii/mu_n, Mob_a = KWr_b*Kii;
5053 const double dMobg_dp = (DKNr_b*Kii/mu_n)*dSe_dp_b, dMobg_dz = (DKNr_b*Kii/mu_n)*dSe_dz_b;
5054 const double dMoba_dp = (DKWr_b*Kii)*dSe_dp_b, dMoba_dz = (DKWr_b*Kii)*dSe_dz_b;
5055 const double gJ = gravity.data()[J];
5056 const double gradSa = -(fsb.
dS_g_dp*grad_u_w_ext_b[J] + fsb.
dS_g_dz*grad_u_n_ext_b[J]);
5057 const double gp_a = grad_u_w_ext_b[J] - rho_a_mass_b*gJ;
5058 const double gp_g = grad_u_w_ext_b[J] + pcp_b*gradSa - rho_g_mass_b*gJ;
5059 ugI -= Mob_g*gp_g; uaI -= Mob_a*gp_a;
5060 const double dgradSa_dp = -(fsb.
d2S_g_dp2 *grad_u_w_ext_b[J] + fsb.
d2S_g_dpdz*grad_u_n_ext_b[J]);
5061 const double dgradSa_dz = -(fsb.
d2S_g_dpdz*grad_u_w_ext_b[J] + fsb.
d2S_g_dz2 *grad_u_n_ext_b[J]);
5062 const double dgpg_dp = dpcp_dp_b*gradSa + pcp_b*dgradSa_dp - drgm_dp_b*gJ;
5063 const double dgpg_dz = dpcp_dz_b*gradSa + pcp_b*dgradSa_dz - drgm_dz_b*gJ;
5064 dugp -= dMobg_dp*gp_g + Mob_g*dgpg_dp;
5065 dugz -= dMobg_dz*gp_g + Mob_g*dgpg_dz;
5066 duap -= dMoba_dp*gp_a + Mob_a*(-dram_dp_b*gJ);
5067 duaz -= dMoba_dz*gp_a + Mob_a*(-dram_dz_b*gJ);
5069 ug_b[I]=ugI; ua_b[I]=uaI;
5070 dug_dp_b[I]=dugp; dug_dz_b[I]=dugz; dua_dp_b[I]=duap; dua_dz_b[I]=duaz;
5075 double F_n_dot_n = 0.0;
5076 for (
int I = 0; I < nSpace; I++) {
5077 F_n_dot_n += (Ag*ug_b[I] + Aa*ua_b[I]) * normal_b[I];
5080 double Sval_p_b = 0.0, Sval_z_b = 0.0;
5081 for (
int I = 0; I < nSpace; I++) {
5082 Sval_p_b += (dAg_dp*ug_b[I] + Ag*dug_dp_b[I] + dAa_dp*ua_b[I] + Aa*dua_dp_b[I]) * normal_b[I];
5083 Sval_z_b += (dAg_dz*ug_b[I] + Ag*dug_dz_b[I] + dAa_dz*ua_b[I] + Aa*dua_dz_b[I]) * normal_b[I];
5089 double Kw_rep = 0.0;
5090 for (
int ii = 0; ii <
nnz; ii++) Kw_rep = fmax(Kw_rep, fabs(KWs_eN[ii]));
5091 const double a_n_scale = (Ag*KNr_b/mu_n + Aa*KWr_b) * Kw_rep;
5092 const double penalty = ebqe_penalty_ext.data()[ebNE_kb] * a_n_scale;
5095 F_n_dot_n += penalty * (u_n_ext_b - bc_u_n_ext_b);
5107 for (
int i = 0; i < nDOF_test_element; i++) {
5108 const double test_i_dS = u_test_trace_ref.data()[
5109 ebN_local_kb * nDOF_test_element + i] * dS_eb;
5111 elementResidual_n_eb[i] += F_n_dot_n * test_i_dS;
5114 for (
int j = 0; j < nDOF_trial_element; j++) {
5115 const double trial_j_b = u_trial_trace_ref.data()[
5116 ebN_local_kb * nDOF_test_element + j];
5119 double Sgrad_p_b = 0.0, Sgrad_z_b = 0.0;
5120 for (
int I = 0; I < nSpace; I++) {
5121 for (
int ii = a_rowptr.data()[I]; ii < a_rowptr.data()[I + 1]; ii++) {
5122 const int J = a_colind.data()[ii];
5123 const double Kii = KWs_eN[ii];
5124 const double Mob_g = KNr_b*Kii/mu_n, Mob_a = KWr_b*Kii;
5125 const double gNjJ = u_grad_trial_trace_b[j * nSpace + J];
5127 const double dFdgp = Ag*(-Mob_g*(1.0 - pcp_b*fsb.
dS_g_dp)) + Aa*(-Mob_a);
5129 const double dFdgz = Ag*(-Mob_g*(-pcp_b*fsb.
dS_g_dz));
5130 Sgrad_p_b += dFdgp * gNjJ * normal_b[I];
5131 Sgrad_z_b += dFdgz * gNjJ * normal_b[I];
5135 double jac_nn = Sval_z_b * trial_j_b + Sgrad_z_b;
5136 double jac_nw = Sval_p_b * trial_j_b + Sgrad_p_b;
5140 jac_nn += penalty * trial_j_b;
5145 elementJacobian_n_n_eb[i][j] += jac_nn * test_i_dS;
5146 elementJacobian_n_w_eb[i][j] += jac_nw * test_i_dS;
5152 for (
int i = 0; i < nDOF_test_element; i++) {
5153 const int eN_i = eN * nDOF_test_element + i;
5154 const int gi = u_l2g.data()[eN_i];
5155 globalResidual.data()[offset_n + stride_n * gi] += elementResidual_n_eb[i];
5156 if (have_gas_budget) {
5157 gas_budget_node.data()[4 * numDOFs_n + gi] += elementResidual_n_eb[i];
5158 gas_budget_node.data()[5 * numDOFs_n + gi] += elementResidual_n_eb[i];
5160 for (
int j = 0; j < nDOF_trial_element; j++) {
5162 + i * nDOF_trial_element + j;
5163 globalJacobian.data()[csrRowIndeces_n_n.data()[eN_i]
5164 + csrColumnOffsets_eb_n_n.data()[ebN_i_j]]
5165 += elementJacobian_n_n_eb[i][j];
5166 globalJacobian.data()[csrRowIndeces_n_w.data()[eN_i]
5167 + csrColumnOffsets_eb_n_w.data()[ebN_i_j]]
5168 += elementJacobian_n_w_eb[i][j];